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1 Introduction 


This final report presents the results of an almost four year effort on finite element modeling 
of the static and dynamic behavior of anisotropic materials with particular emphasis on 
single crystal alloys used in the manufacturing of high-performance turbine blades. This 
effort was motivated by a lack of finite element software capable of representing stress and 
strain as second order tensors. 

During the course of the project, various formulations for two— and three-dimensional 
hybrid finite elements were developed and implemented into the SPAR finite element code. 
These formulations were tested and compared with displacement-based approaches for both 
static and dynamic problems. 

In an extension of the original statement of work, a sensitivity analysis of experimental 
results for anisotropic materials subject to misalignments and other errors was conducted. 
As a results of this study, a formulation, numerical procedure, and computer software were 
developed for the calculation of material constants for anisotropic materials and, moreover, 
for the optimization of strain gauge and material axis orientations of tensile test specimens 
for highly anisotropic materials. Moreover, as an additional task, a similar procedure and 
software were developed for the evaluation of stresses by means of strain measurements and 
for the optimization of the orientation of strain gauges in this case. 

This report presents a summary of the technical effort over the course of this project. 
The report is divided into several sections: In Section 2, general hybrid finite element for- 
mulations are studied. Then, in Section 3, two-dimensional finite elements based on hybrid 
formulations are developed and tested numerically. In the next section, a dynamic analysis, 
using hybrid finite elements is discussed. This analysis is followed by the formulation of 
three-dimensional hybrid finite elements (Section 5). In Section 6, alternate hybrid stress el- 
ements are studied and verified numerically. As a conclusion of the theoretical developments 
presented, several numerical examples cases for the SSME turbine blades are presented in 
Section 7. This section is followed by the formulation of numerical procedures for the op- 
timization of the orientation of material axes and strain gauges in experimental tests for 
anisotropic materials (Section 8). Several numerical examples are presented to illustrate the 
performance of the procedure. 

In the appendices, updates to the SPAR code manual are presented. In four separate 
volumes, theory and user’s manuals for the two experiment optimization codes — OPTAM— C 
and OPTAM-S — are included. 
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General Hybrid Element Formulation 


2.1 Introduction 

This project began with the study of a variety of methods for modeling vibrations of 
anisotropic elastic blades with particular emphasis given to hybrid finite element formulations 
and the feasibility of using helical shell elements. Apparently, some authors have attempted 
to model isotropic curvilinear turbine blades using shell elements of variable thickness based 
on a shell theory for surfaces generated along a helix. One advantage in this approach is 
that models with relatively lew degrees of freedom can yield quite acceptable results. A 
disadvantage, of course, is the limited flexibility inherent in shell models for capturing the 
effects of boundary conditions, three-dimensional stress states, etc. 

Several alternative hybrid element formulations were investigated. Purely qualitative 
studies were made, the objective being to assess a priori properties of various elements with 
regard to 

1. accuracy in computing stresses 

2. ease in handling anisotropic properties 

3. numerical stability in the presence of strong anisotropy 

4. compatibility with the SPAR code structure 

The basic formulation which generated interest is that which yields a hybrid element 
from assumed displacement fields. Starting with the principle of minimum potential energy 
with displacement continuity conditions as constraints, the boundary tractions then appear 
as Lagrange multipliers. The resulting functional is of the form 

v — X] [J v (^)Cijki£ij{u)eki(u) - F,vrj dV - Tiu,ds + T t u,ds^ (2.1) 

Here standard notation is used: V e is a typical element. C,jki are the elastic constants, 
Sij the strains, F, the body forces, T x the prescribed tractions on dV\ e and u x the prescribed 
displacements on dV 2e , where the T x in this integral are the Lagrange multipliers. 

Approximations take the form 


u = A/3, T = 4>q 

where 0 and q are vectors of undetermined parameters. The strains are then 

£ = B8 


( 2 . 2 ) 


(2.3) 


2 


Denoting 


Hx= [ B T CBdV H 2 = [ A T FdV 

Jv c JV' 

H 3 = [ A T 4>dV H 4 = [ u T <f>dV 

JdV 1. JdV 7, 


(2.4) 


a minimization of i r yields 
Hence, 


(2.5) 

( 2 . 6 ) 


W 


/here 


P=h;\h„ + h 2 ) 

“ e 

F. = Hj(H^') T H 3 
V, = Hj irjH;, + //, 


Equilibrium and continuity of interelement displacements is achieved in the discrete model 
whenever 

F e q = U e (2.S) 


The hybrid element stiffness matrix is then 

K - F7 l 


(2.9) 


Hybrid element formulation 

Six stress and displacement type elements were chosen for studies of stress accuracy and 
utility in vibrational analysis. Two separate finite element formulations were chosen. 

1. Hellinger-Reissner Formulation 

The functional which assumes a stationary value is given by: 


7T r 



1 

9 


<r T Sa + <j t {Du) 



(2.10) 


Here a is the vector of stress components, D is the differential operator defining the strain- 
displacement relations (s = Du), with u the displacement vector, and T the boundary 
traction. 

The -displacements u are not assumed to be compatible and the equilibrium equations 
are not satisfied identically but are brought in as constraints, so that the elements are more 
“flexible” and the number of stress terms required to suppress zero-energy modes is reduced. 
The elements also becomes less sensitive to changes in the reference coordinates [12]. 
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The strain-displacement relation is given by 

£ = Du 


( 2 . 11 ) 


where the displacement u consists of a compatible part u c and an incompatible part u a , 
which could be a bubble function vanishing on the boundary. 

Now, 

J v a T (Du a )dV = - J^(D T a) T u a dV + ^{N a) T ujs ( 2 . 12 ) 

where Na is the trace of a on the boundary, N being a matrix of direction cosines of the 
unit exterior normal to dV . Thus, if we set u a = u — u and Ncr = 7\ we have, 

*, = / [-Vs <7 + t t {Du c ) - (D t <j) t u\ iV (2.13) 

Jv 12 J 

We see that the equilibrium equations 

D T a = 0 (2-14) 

appear as a constraint with the u a being the Lagrange multiplier. 

In the finite element implementation, we assume 

a = P3 (2-15) 


where 

and Ci$ are row vectors. 
Also, let 

and 

from which 
and 
so that 


C2 

u c = Nq 
u a = LX 

Du c = Bq (B = DN ) 

D T a = E0 (E = D t P) 

?r r = -1 fHP + pGq - 0 T Rx A 



(2.16) 

(2.17) 

(2.18) 

(2.19) 

( 2 . 20 ) 
( 2 . 21 ) 
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where 


and 


H = / P T SPdV 
Jv 

(2.22) 

G= [ P T BdV 
Jv 


Rt = / E T LdV 

Jv 

(2.23) 

with respect to /? and A, w'e get 


0 = H~ l (Gq-R rA) 

(2.24) 

Rfp = 0 

(2.25) 


and 


Now, the strain energy as expressed in terms ot a stiffness matrix A and the geneialized 
coordinates q is given by 


U = \q T I<q 

f 1 


U = 


Jv 2 


) 


a 1 edV 


- I <j T S adV 
2 Jv 


(2.26) 


U = -fHd 


Eliminating A from equations (2.24) and (2.25), substituting 3 into the expression in 
equation (2.26) which becomes 

(2.27) 

(2.28) 
(2.29) 


gives 


where 


K = G t MG - G t M R x (Ri M Ri)~ 1 Ri M G 
M = H~ t 


The inversion of H becomes easier if the same C t ’s are used for all of the normal stress 
components. For example, in the case of a three-dimensional isotropic solid, when the stress 


5 


terms are not coupled, 


H = 1/E 


<t> i ~v<t> 1 

- V<t>\ <t>\ ~ v( t> 1 

- i/d>i -isd>i 4> i 


2(l + i/)^ 


2(1 + l /) 4> 5 


2(1 + u)d> 6 


(2.30) 


where 


&= / ^ 

Jv 


dV 


(2.3i ; 


77 can be easily inverted. 

2. Hu- Washizu Formulation 

Another approach to the development of mixed elements is to use the extended variational 
principle of Hu and Washizu. The Hu- Washizu variational functional is given by. 

’I 


xhw — 


IV L 


-2 £ T Cc - cr T z + cr T {Du) 


dV 


■fdV 


T 1 ( u - u)dS 


(2.32) 


where 


C = S~ l (2-33) 

and the independent variables are the stresses <r, the strains S , element displacements u. and 
boundary displacements u. 

In the finite element formulation, both <j and £ are approximated by the same shape 


function: 

sex 

II 

b 

(2.34) 

and 

£ = Pot 

(2.35) 

The strain energy is then 

U = f -o T edV = ll3 T Aa 

(2.36) 

where 

h 2 2 



A = 


B x 

B 2 


(2.37) 
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and 


Bi — J Pf P t dV (2.38) 

The best choice for the reference coordinates should be such that the B{ s are diagonal 
matrices and the inversion of H becomes easy. 

Using the same approximation for u as in equations (2.17) and (2.18) and following a 
similar line of logic results in the following: 

khw = -a T Ja — 3 t Aa + fi T Gq — d T R\X (2.39) 

where 

J= / P T CPdV (2.40) 

Jv 

Setting the first variation of z hw with respect to ft and a to zero gives 

Aa = Gq-R x A (2.41) 

and 

Ja = jf A (2.42) 

Substituting the expression for a from equation (2.41) and (2.42) into the strain energy 
solution (2.36), 

K = G t MG - G t MR x (R fMR^-'RjM G (2.43) 

where 

M = A~ x J A~ l (2.44) 

Although the finite element formulation using the Hu-Washizu functional exists, all ap- 
plications are derived with regards to shell/plate theory rather than analysis of continua. 

Pian [15] has worked out solutions for three-dimensional brick elements using the mod- 
ified Reissner principle. He has analyzed the 8 node hexahedral solid and the 20 node 
hexahedral solid. 


Hybrid Element Shortcomings 


Stress-hybrid elements possess some shortcomings which must be overcome in order for them 
to be useful. The first major shortcoming of hybrid stress elements is that they can possess 
zero energy modes which can have debilitating effects on eigenvalue problems. Babiiska, 
Oden, and Lee (CMAME, 1978) have shown that such methods are stable and convergent 
only if a global LBB-condition is statisfied, i.e., conditions of the type 


“Ill-Mil < su p 

u 



(2.45) 
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where A is a multiplier on the constraint of continuity of displacements u . across (mt «ejement) 
boundaries and ||| • |||, || • 111 are appropriate norms (for details, set i ])■ h f 

mesh refinement. 


To overcome this problem, the method must possess an LBB-parameter 

0 < a = constant independent of h 


(2.46) 


4 nuick check for a necessary condition for stability is to calculate the rank of the matrix 
associated with the constraint term (e.g., § A uds). For the Bellinger- eissnei oimu a »on 
““ energy modes may normally be suppressed when the number of 
to or larger than the total degrees of freedom minus the number o rig 
freedom. However, these stress terms must be chosen carefully because they do 
contribute to the suppression of zero energy modes. 

For the Hu-Washizu formulation, in the eight node solid, the minimum number of s ress 
terms required to eliminate zero-energy modes is (S x 3 - 6 = IS). Using these terms, ,.e„ 

<j x = 8\ “1" fay "t" 4* 


G 7 = 


(2.47) 


T xy = /3i3 + P 


14 * 


Ty Z 4” 1^15 "4 8\6 X 


T xz = /?i7 + PisV 


the bending of a beam was analyzed, and the results were not good. 

If, however, the stresses were to be decoupled as suggested, additional terms would be 
required to eliminate the zero-energy modes, in fact, 21 0 s would be necessary. 

By using three internal displacement parameters (A’s), the number of 0's may be reduced 
to 18. Results obtained with this formulation were quite good. 

For two-dimensional problems, especially plane stress, Spilker [24] has shown th ® 
of an eight noded quadrilateral with complete stress terms. His paper also considers 4 noded 
quadHlaterah^using 5 0's for the stress terms and a bilinear 
Although the results were good,' the element was sensitive to change 

dinates. . , , 

The second major shortcomrng for hybrid stress elements is that these eU ™ents ma, 
stiffness matrices that are not invariant under a change ,n local coord, nate systems. Indeed, 
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the original Pian-hybrid-stress element (now in SPAR) can yield severely different stresses 
when simple changes in the local coordinate system occur. 

Invariance of the stiffness matrices may be guaranteed by introducing appropriate bubble 
functions” which serve as multipliers on the equilibrium constraint, 

D T a = 0 ( 2 - 4S ) 


Completeness in all modes of the polynomial expansion is also important, to ensure 
invariance with respect to the chosen reference coordinates. Another method of achieving 
invariance is to use local coordinates. 

In the modified Hellinger-Reissner formulation, the equilibrium equations are introduced 
as a constraint and do not have to be satisfied pointwise. The degree of satisfaction depends 
on the number of internal displacement parameters A that are used. 

The third and one of the most serious shortcomings of stress-hybrid elements is that, in 
their present form, they are incapable of yielding a consistent approximation of the kinetic 
energy in an element. This is a deficiency often “swept under the rug in discussions of 
applications of hybrid elements to problems of structural vibrations. Traditionally, mass 
matrices are calculated using displacement or velocity approximations which are comp ete y 
independent of those used for (or resulting from) calculations of the element stiffness matri- 
ces. To overcome this deficiency, it has been observed that the kinetic energy in a body . 
can be written , , r , 

(2.49) 


b 




where p is the mass density and p is the momentum. The momentum is related to the stress 
by p = D r a (2-50) 


Hence, it may be possible to define a “consistent approximation” of the kinetic energy which 
should produce more accurate approximations of mode shapes and frequencies. 
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3 Two-Dimensional Element Definition 


3.1 Definition of Element Matrices 


Using the Hellinger Reissner formulation and introducing the equilibrium conditions as a 
constraint into the complementary energy functional and modifying produces the following 
functional: 

= <r T SadV-[ a T (Du)dV + f u T Tds\ (3.1) 

„ 12 JVn JVn JdVn J 

where cr,u are the stress and displacement vectors, 5 is the compliance matrix, D is the 
differential operator, V n is the volume of the nth element, and T is the prescribed stress on 
the boundary S an . 

The stresses are interpolated in terms of the stress parameters 3 and the polynomial 
shape functions P, i.e., 

<j = P/3 (3-2) 


so that equilibrium is satisfied, i.e., 


For isoparametric elements, 


x = 


V = 


D T a = 0 


i 

T. 

i 


(3.3) 


(3.4) 


. = y 

i 


where (£,f?,0 form the parent plane and iV,(£, ??, s) are the appropriate shape functions. 
The displacements are interpolated in terms of the shape functions as 


u = N(£,Ti,q)q 


(3.5) 


where q are the element nodal displacements. 
Then the strain is given by 


s = Du = B(Z,Tj,<;)q = 


where B = DN(x,y,z ) and |J| is the Jacobian of the transformation. 


(3.6) 
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Substituting equations (3.2), (3.5), and (3.6) into (3.1) and defining 


- - IXj 

f P T SP\J\did V dq 
-1 


G - IXj 

' P T B’d^dqdq 
-1 

( 3 . 7 ; 

Q = j\ ^ N T Tds 


and noting that 

dV = 

: \J\d^drjdq 

(3.8) 

gives 

n ^ w 

'H3 - 3 T Gq + C/} 

(3.9) 

Equating the first variation of i r mc to zero i 

in each element results in 


id = H~ l GLq* --- H~ l Gq 

(3.10) 


where q* are the global displacements which are related to the element displacements by the 
Boolean matrix L a . 

By substituting (3.10) into the expression for b7r mc , a new expression for 67r mc may be 
written as 

fame = 6q mT L T K Lq* ~J2 lT q\ = ° (3.11) 

In n ) 

where Q is the element load vector and K is the element stiffness matrix 

K = G t H~ 1 G (3.12) 

The equation (3.7) can be numerically integrated, but care should be taken to ensure 
that the proper order of integration is used. 

It has been proved in the paper by Spilker that if the assumed stresses are complete 
polynomials, the element stiffness will be invariant to general rotation and translation. To 
reduce the numbers of /Ts, the equations of equilibrium are applied to the assumed stress 
terms. However, lack of invariance does not imply poor element performance. In isopara- 
metric elements, however, as the choice of the local orthogonal system is not unique, the 
element must be inherently invariant for good results. 

As expected in plane elements, the optimal sampling points are the Gaussian points of 
integration for the stresses. However, the results are still better than those in the assumed 
displacement method, as the stresses in the hybrid method satisfy equilibrium conditions 
which relate the stress gradients. 


11 


3.2 Examination of Different Element Models 

The four-node linear displacement model (see Fig. 3.1). 

The shape functions iV,(£, if) for the intra-element displacement field correspond to those 
for four-node bilinear functions. 

Using complete polynomials for the stress approximation we have for a linear field 

a x = Pi + $2 y + 06 x 

a y = fa + fax + fay ( 3 - 13 ) 


a xy = fa- fax - fay 


which is variant. 

The minimum number of fas 
rank is 


required, yp, as a necessary condition for correct stiffness 


yp > Vd.o.[ - Vr.d.oS 


(3.14) 


where 


= number of element d.o.f. 

= number of rigid body d.o.f. 


(3.15) 


Spilker, et al. compares the performance of two elements, one with 7 fas (PH7L) and one 
with five fas (the minimum required by equation (3.14) - PH5L). 

The optimal sampling point for both these elements is the centroid of the element. The 
dimensions of the element are a and b. The stiffness matrix for the element PHiL corre- 
sponding to a set of generalized displacement parameters a, (for plane stress) is 


1 



0 


V 

(W^) 


0 


0 


0 \b 2 (l + h) 0 


0 0 


K a = 4 Eab 


v 

(T^) 


o 


0 ' o 

0 0 


(l-*' 2 ) 


0 


0 

0 


-a 2 (l + fi) 0 
1 1 

° 2(1 + fa . 


(3.16) 
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Figure 3.1: The four node linear displacement model. 
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i.e., the displacements are interpolated as simple bilinear polynomials of x,y in terms of the 
eight displacement parameters. When the strain displacement relations (3.6) are used, to 
relate £ to a, only five a’s remain as the constant terms fall out. These five a’s correspond 
to stretching and bending along the two axes and pure shear: 

Qi - |a(-ui +u 2 + u 3 - w 4 ) 
a 2 = \ab(u x -u 2 + u 3 — ti 4 ) 


a 3 = \b(-vi - v 2 + v 3 + u 4 ) 
a A = \ab(vi -v 2 + v 3 - v 4 ) 

c*s = \db [6( — Ui — u 2 + U3 + u. 4 ) + a( — ui + v 2 + t’3 — U4)] 

where all of the displacements correspond to the local coordinate system used. 


The terms f\ and f 2 are 

, = [(a/6) 2 -t,] 2 

Jl [1 - i/ 2 + (a/6) 2 2(l + 1/)] 


[jb/al ~ H 2 

[1 - v 2 + (b/a) 22 (l + v)} 


(3.17) 


(3. IS) 


K a is definite and hence, PH7L has no spurious zero energy modes. 

In PH5L, fi and f 2 are different and, therefore, cause a difference in the analysis of the 
bending problem. Element PHhL shows no spurious zero energy modes at 9 = 0. But since 
K depends on 9 the stiffness rank should be explored. 

Computing the G matrix at an arbitrary angle 9 , in terms of the general displacement 
parameters 


where 


4ab 


0 Uab 3 — a 3 bc 2 ) 0 — |a6 3 C! 0 


G° = 


4a6 


0 |a 3 6ci 

0 0 


0 ~ ab 3 c 2 ) 0 

0 0 4a6 


sc(c 2 — s 2 ) s 2 c 2 

s 6 + c 6 ’ ° 2 c 6 + s 6 


(3.19) 


(3.20) 


14 


with s = sin 9 and c — cos 9. 

It is fend that the matrix G becomes singular for J = S 

modes, so that PH5L is not completely invariant. Thus, this elemen 

grea <\ pure bending problem was solved using each of these two elements; the results are 

j rJ 3 0 PH7L leads to a “stiff” solution that converges only when more 
summarized m rig- . • , i:i.» prr^T and gives good results at 

than ten elements are used. However, it is invariant, unhke PH5L. and gives 0 

all 0’s. 

The eight node quadratic displacement model (see Fig. 3.3) 

The shape functions here correspond to those for eight node isoparametric serendipity ele- 
ments (for the displacement approximation). . „ . 

For this 16 degree of freedom element, a minimum of 13 0's is required. ie sma es 
stress field that satisfies this requirement in equilibrium and invariance, is a cub . 

After equilibrium conditions are satisfied, an 18 0 field is derived 

a x = 0 j + ;3 6 x + 02V + fay 2 + 2 3c>xy + 3 10 x 2 + 3 13 x 3 
+ 3.^14 x 2 y + 3,0i oxy 2 + 3 17 y 3 


) xy 


02 + 04 x + fay T 2 0 n xy + 3 u x 2 + faoy~ + 3fai x V 
+ 0xa y 3 + 3 8 l6 x 2 y + 3 ls x 3 

0 S - fax - 0 6 y - 0n* 2 - fay 2 ~ 2 fa° xy - 3 fa* x ~ y 


( 3 . 21 ) 


,.3 


— 30i 4 xy 2 — Pisy 3 — fa$ x 

which is invariant and has no spurious zero-energy modes. 

To reduce the number of »' s, the Beltrami-Mitchell compatibly conditions (for p ane 
strain) are employed, where (3.22) 

V 2 (cr r + a y ) = 0 V 
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Then, after re-numbering 


°* = 0i + 0ex + 0 2 y + 0 s y 2 + 2 0 9 xy + R 10 x 2 + 0 n x 3 

+ 0iz{3x 2 y - 2y 3 ) + W u xy 2 - 0 l5 y 3 
a y = 03 + 04 X + 0 7 y + 2 0 n xy - 0 8 x 2 + 0 lo (y 2 - 2x 2 ) 

(3.23) 

+ /?i2(3 xy 2 — 2x 3 ) + 0i 3 y 3 + 30 i5 x 2 y — d^x 3 
°xy = 05 ~ 07X - 0 6 y - 0 n x 2 - 0 9 y 2 - 2 Pioxy - 30 u x 2 y 

- 3013 xy 2 - 0 l4 y 3 - 0 15 x 3 

This element (PHloQ) is still a complete cubic with only 15 0's and is invariant and has 
no spurious energy modes. The optimal sampling points for both the PHI 50 and PH180 
elements are the 2x2 Gaussian points. 

Another reportedly good element is the 14 3 element developed and tested bv Hershell. 
However, it is not invariant. 

Some example problems were solved and the results, as well as the conclusions are 
summarized in Fig. 3.4. 

PIT h^ P 1 revi ° US . Sraph gives the results for a P lane str ess problem. The performance of the 

!i v!uio^ nt 13 eXCeI I ent ’ 6Ven When the order of integration is reduced to 3 x 3. However, 
the PHlbQ element stiffness becomes singular when the integration order is reduced. 


3.3 Numerical Experiments 

Some numerical experiments were conducted to compare hybrid stress elements with tra- 
ltional displacement method elements. Results show remarkable deficiencies in traditional 
displacement methods for anisotropic materials. The analysis of a cantilever beam under 
he effect on an end shear load was selected as the test problem, as shown in Fig. 3.5. The 
data for the problem include the following: 
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Figure 3.5: Anisotropic cantilever beam with an end load 
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= Young’s modulus in direction 1 
E 2 = Young’s modulus in direction 2 


i/ 2 = Poisson’s ratio 


G 12 = Shear modulus (independent) 

9 = angle between material and global axes 

For the isotropic case, classical beam theory gives 
deflection: 


We shall use 

W = 1501b., / = 10”, E = 3.10' lb/in 


_ Wl 3 
Utip " 3EI' 


Then 


^tip 


. 02 ” 


the following expression for the tip 

(3/24) 

2 , I = h 3 /12 = \ (3.25) 

(3.26) 


The assumption of plane stress was made. 

(a) Isotropic case 

L Linear quadrilateral elements (4 nodes, Id, meshes shown in Fig. 3.6) 


No. of nodes 

(elements) 

(mesh) 

«ti P (hybrid) /u an ai. 

“tip (disp.)/ “anal. 

( Ur )max (hybrid) 

(^i)max (anal.) 

( )max (disp. ) 

( )max (anal.) 

22 (10) (A) 

0.8815 

0.6789 

0.7837 

0.7036 

33 (20) (B) 

0.8297 

0.7110 

0.7442 

0.7152 

44 (30) (C) 

0.8128 

1.1578 

0.7418 

1.2239 

63 (40) (D) 

0.9690 

0.8950 

0.9064 

0.9176 


2. Quadratic quadrilateral elements (8 nodes, 15/3) 


No. of nodes 

“tip (hybrid)/u a nal. 

^tip (disp.)/ n ana j > 

(“x)max (hybrid) 

(o-x)max(disp.) 

(elements) 

(^x)max (anal.) 

( & x )max (anal.) 

28 (5) 

1.0070 

0.9885 

0.9841 

0.9693 
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Figure 3.6: Finite element meshes used for testing the two-dimensional hybrid stress models. 
A. 10 elements B. 20 elements 

C. 30 elements D. 40 elements 
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(b) Anistropic case 

1(a) Linear quadrilateral elements 

^1 - 3 x 10 7 psi, E 2 = 3 x 50 5 psi,v 12 = 0.3 
G 12 = 1.1538 x 10'psi, 9 = 30° 


No. of nodes 
(elements) 

| “tip (hybrid) 

“tip (displ.) 

( Ox) m (hybrid) 

(^r)m (displ.) 

63 (40) 

0.1136” 

0.6271” 

6994.3 psi 

2353.4 psi 


1(b) Quadratic quadrilateral elements 

No. of nodes = 28,9 = 30°, E l /E 2 = 100 


“tip (hybrid) 

| «ti P (displ.) 

(^x)m (hybrid) 

1 {vx) m (displ.) 

r xy (hybrid) j 

r xy (displ.) 

0.1411” 

0.S799 

6827.0 psi 

6452.8 psi 

697.8 psi 

419.2 psi 


Similarly, lor the same values ol E\, E 2 , G'i 2 , z/ 12 , but for 9 = 45° and 9 = 0°, the above 
cases were run. The results obtained were: 

2(a) Linear quadrilateral elements 6 = 45° 


“tip (hybrid) 

“tip (displ.) 

(<7r)m (hybrid) 

(a x )m (displ.) 

Try (hybrid) 

T xy (displ.) 

0.4049” 

0.4131” 

6867.2 psi 

10174 psi 

1362.8 psi 

733.3 psi 

2(b)Quadratic quadrilateral elements 6 - 

= 45° 



“tip (hybrid) 

“tip (displ.) 1 

(&x )m (hybrid) 

(o'x)m (displ.) 

r xy (hybrid) 

r xy (displ.) 

0.4964” 

0.8831” j 

6682.0 psi 

7844.6 psi 

876.3 psi 

953.2 psi 
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3(a)Linear elements 0 = 0° 


«tip (hybrid) 

u tip (displ.) 

(Vx)m (hybrid) 

{°x)m (displ.) 

T xy (hybrid) 

r xy (displ.) 

0.02009” 

.018359” 

6859.3.0 psi 

6314.5 psi 

214.6 psi 

1019.4 psi 

3(b) Quadratic elements 0 = 0° 




«t. P (hybrid) 

u tip (displ.) 

{o x )m (hybrid) 

(<7 r ) m (displ.) 

r xy (hybrid) 

T xy (displ.) 

0.02022” 

0.02007” 

6804.1 psi 

6687 psi 

215.2 psi 

248.2 psi 


The analytic solution for the case with 0 = 0° is 

<7 X = 6971.4 psi 

r xy = 225 psi 

So. the errors in the hybrid element model are: 

<7*: 2.4% ; r xy AA% 

and in the displacement model, the errors are 

cTj.: 4.2%(at best ) ; r rv :ll.l% (at best) . 

Next, the case of a tapered cantilever beam was considered (see Fig. 3.7). Only the 
elements that gave reasonably accurate results were used (the quadratic quadrilateral hybrid 
stress elements), which were compared with the regular displacement type elements. The 
tip displacement, as determined by elementary beam theory for a uniform load IF /length is 
u = 4.08 x 10 3 inches. Beam theory gives for the bending stress at any point is a x = h—i 
at that section, and is maximum at the supported end. 

(a) Isotropic case - using quadratic quadrilateral elements (8 nodes, 15 /3), 


No. of nodes 
(elements) 

Utip(hybrid) 

^tip(anal.) 

«ti P (disp.) 

a^hybrid 
a x anal. 

<7 r disp. 
a x anal. 

45 (10) 

1.0298 

1.0164 

1.0180 

1.0323 
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(b) Anisotropic case - 8 nodes, 15 (3 elements 


E\— 3 x 10'psi ; E 2 /Ex = 1/100 


No. of Nodes 
(elements) 

u tip (hybrid) 

^tip (displ.) 

a x hybrid 

<7 r displ. 

9 = 45° 
45 (10) 

0.1060” 

0.2012” 

2068.0 psi 

2636.8 psi 

9 = 30° 
45 (10) 

0.0315” 

0.03876” 

1968.5 psi 

2076.0 psi 

9 = 0° 
45 (10) 

0.00423” 

0.00419” 

2027.3 psi 

1999.0 psi 


Discussion of Numerical Experiments 

For isotiopic materials, when the analysis was conducted using four-node quadrilateral ele- 
ments, neither the hybrid nor the displacement models captured the cubic variation of the 
displacement until the mesh was sufficiently refined. However, the stresses from the hybrid 
model were better, though only marginally, since only a linear variation in the stress was 
used. When quadratic elements were used, the hybrid element model gave excellent results, 
both tor stresses and displacements, as the stress shape function was a complete cubic and 
was thus able to approximate the bending moment and shear stress very accurately. 

For anisotropic materials, although the results obtained using linear and quadratic hy- 
brid models were close, the displacement models gave significantly different results, in both 
displacements and stresses. This behavior is especially significant when the material axes 
are inclined to the global axes. While the hybrid element model continues to give excellent 
sti esses, the displacement model gives very poor results. 

For a tapered cantilever beam, again, the hybrid model gave better displacements than 
the displacement model, especially for anisotropic materials with material axes that do not 
coincide with the global axes. 
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4 Vibrational Analysis 

4.! A Variational Principle for Dynamic Analysis 


The equations 


for dynamic equilibrium are given by. 


ffijj + Fi = r ‘ 


m V at all t 


where n = inertia impulse vector. 
Equation (4.1) may be written as 


[/w + ^ - r 'H = 0 


Hence, another way 


of writin 


cr this equation would be 


j + /. _ r * 


= ffjj and /. = Fi. tCTSor T ‘: 


IS 


referred to as the impulse tensor 


where Tq — ana * «• - ' J surface impulses. 

The boundary conditions are prescnbed veloat.es and surtace P 

mutually exclusive regions, i.e., 


U| 


on s. 


and 


t.jUj = bi on s 1 


where n, is the component ot the unit normal vector on s,. 

The strain-displacement relation is 

£ij = 


and the kinetic energy function is 


U 


= [ u\( r <) dV 

Jv 


while the complementary strain energy is 

U‘ = I u\(r t] )dV 

Jv 


where 


In linear elasticity, 


du‘/dT ZJ = £ij 




(4.1) 

(4.2) 

(4.3) 

field. 

specified on 

(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 

(4.9) 

( 4 . 10 ) 
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Corresponding to Hamilton’s principle, we can define a functional 


as 


7T C 


- f U( r i) - - f tiUids 

L J 3u 


dt 


(4.11) 


where the assumed stresses must satisfy the equations of motion and the traction (impulse) 
boundary conditions of s T . 

In the finite element formulation, equilibrium must be satisfied in each element and on 
inter-element boundaries so that 


7T = 


-/ E Up - Up( a ij,}) - [ tiUids - f u t t,ds 

-n ^ JS11P Jskjp 


dt 


( 4 . 12 ) 


where s u p denotes the part of the surface of Vp on which velocities are prescribed and syp 
denotes the interelement boundary of Vp, the particular element 


S P — s uP + S T p + SjVP 


so that 


dt 


(4.13) 

(4.14) 


T = / E Up - u; - ! tiUids - / udids 

•hi L J Sp J S r p 

The admissible surface velocities and the impulses must satisfy the following conditions: 

1. Prescribed at arbitrary times t\ and i 2 * 

2. Continuous first derivatives for t X j , continuous m. 

3. Equilibrium conditions (dynamic/boundary conditions). 

4.2 Formulation of Element Matrices for Dynamic Analysis 

For dynamic analysis using the hybrid/displacement models, a consistent mass matrix was 
generated with a kinetic energy term that was introduced into the energy functional, so that 

Kdyn = - J v ° T SsigdV - j^sig 7 DudV + u T Tds - J pii T udV ( 4 . 15 ) 

When discretized, this becomes 

^dyn = E J v ScdV — <r T Dudv + J u T Tds — j pu T udV ( 4 . 16 ) 

If u = Nq, the kinetic energy term becomes 

J v pq T N T NqdV = q T Mq ( 4 . 17 ) 
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where 


M = j v N T P N\J\did v (4.18) 

is the mass matrix for each element. 

After the element mass matrices are assembled to form the global mass matrix, the 
generalized eigenvalue problem 

Kx — A Mx (4.19) 

is solved for the first few eigenpairs. 


4.3 Numerical Experiments for Hybrid Stress Element Vibra- 
tional Analysis 

Problem Definition 


A cantilever beam is analyzed tor its first tew natural frequencies (eigenvalues) and mode 
shapes using both the assumed-displacement and the hybrid— stress method. 

A consistent mass is generated and the generalized eigenvalue problem 


Mq + Kq = 0 


(4.20) 


is solved for its eigenpairs, which are the natural frequencies and mode shapes of the physical 
system. Since the size oi the matrices is not very large, a solver from IMSL that determines 
the eigenvalues and eigenvectors is used instead of the sub-space iteration scheme suggested 
by Bathe [2], 


For Bernoulli-Euler beams composed of isotropic materials (neglecting the effect of shear 
deformation and rotatory inertia since only the first two modes will be considered where 
the correction introduced as a result of these effects is small), the equation of motion of 
transverse vibration is 


<£_ 

dx 2 



a d 2u 


(4.21) 


where 

u = u(x,t) is the transverse displacement 


A — area of cross section of the beam 


x = axial distance from the point of support 

p = mass density of the material 

/ = centroidal moment of inertia of the cross section 

The boundary conditions for the cantilever are: 

At the fixed end: 
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(4.22) 


At the free end: 


du 

u = 0 and — = 0 
dx 


d 2 u d 3 u 

d?=° and i?=° 


(4.23) 


Substituting the boundary conditions into the general solution, we get three homogeneous 
linear algebraic equations which would give a non-trivial solution only if the determinant ol 
the coefficients vanishes, are derived 


1 + cos XL cosh XL + 1 = 0 


(4.24) 


which is the characteristic equation whose roots are the eigenvalues X r times length L. A 
numerical solution alone exists for the above equations, determined by Craig and Bampton. 
The first few values are 


X 1 L = 1.8751 
X 2 L = 4.6941 

and the natural frequencies for the cantilever are given by 

(A rJ L) 2 CE/N* 


so 


that 


L 2 \pA) 


3.516 / ElV 




L 2 


- 


,P A ) 


and 





22.03 (Eiy 

U~ \Ja) 


(4.25) 


Substituting the numerical values for the given problem results in the following: 

wi = 17.58Hz , to 2 = 110.15Hz 

The mode shapes are given by 

V T {x) = cosh(A 4 x) — cos(A r x) — A: r [sinh(A r x) — sin(A r x)] 

where 

' cosh(A 7 .L) + cos(A r L) 
sinh(A r T) + sin(A r T) 

as in Craig. 

The first two mode shapes for a cantilever in free vibration are shown in Fig. 4.1. 


(4.26) 

(4.27) 

(4.28) 

(4.29) 

(4.30) 

(4.31) 
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Linear Element Results 


The meshes that were used for the static problem are used here, with the number of elements 
varying from 10 to 80. 

The normalized natural frequencies (uq anal./uq f.e. and u 2 anal ,/u 2 f.e.) are plotted 
against the number of elements, and are shown in Figs. 4.2 and 4.3. 

For the isotropic case, the hybrid model converges to the analytical solution faster than 
the assumed displacement model. The mode shapes however do not seem to vary much, as 
seen in Fig. 4.4 (for the mesh with 80 elements). 

The material model chosen for the anisotropic case is that of cubic syngony with the same 
properties as in the static case, and the 40 element mesh. The first two natural frequencies, 
for various angles of rotation of the material axes, using both finite element approximations 
are tabulated in Table 4.1 


^lamcrxT 

(3 7> 



£J ro ?C 


O . 2 tizrr^rj^b 

Figure 4.2: Normalized natural frequencies (mode 1 analytical solution vs. mode 1 finite 
element solution) for various linear element mesh sizes. 



1C 2 O . 30 <o 50 6o 70 9o 

4 °0* cj ^J^rrLerOb 


Figure 4.3: Normalized natural frequencies (mode 2 analytical solution vs. mode 2 finite 
element solution) for various linear element mesh sizes. 
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Figure 4.4: Linear quadrilateral finite element meshes used for the numerical experiments 
with hybrid stress elements and vibration analysis. 
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Table 4.1 


Natural frequencies for an anisotropic cantilever beam 
using 40 linear elements for different material axes’ orientations 


Orientation of Axes (hybrid) 


0 

0 

14.351 

Hz 

00 

0 

0 

13.307 

Hz 

45° 

12.873 

Hz 

60° 

13.307 

Hz 

<0 

O 

0 

14.351 

Hz 


uii (displ.) 

u 2 (hybrid) 

(displ.) 

14.785 Hz 

85.950 Hz 

88.341 Hz 

14.726 Hz 

80.673 Hz 

88.0945 Hz 

14.058 Hz 

78.379 Hz 

85.011 Hz 

14.547 Hz 

80.673 Hz 

87.318 Hz 

14.247 Hz 

S5.950 Hz 

S5.934 Hz 


From the above table it is observed that the hybrid model gives identical results for a 
lotation of 90 and no rotation of the material axes, and for 30° and 60° rotations of the 
axes. The displacement method however gives results that vary, even though the moduli E\ 
and E 2 are equal. 

Quadratic Element Results 

The natural frequencies and mode shapes of the isotropic and anisotropic cantilever beams 
are now calculated using an eight noded finite element mesh with the number of elements 
varying from 3 to 20. The meshes used are the same as those for the static case. 

The normalized natural frequencies {u>i anal./uq f.e. and anal./u.^ f.e.) are plotted 

against the number of elements, and are shown in Figs. 4.5 and 4.6. Again, it is observed that 
the hybrid model converges to the analytical solution faster than the assumed-displacement 
method. The mode shapes however are very similar in both models, except for the maximum 
amplitude (when 20 quadratic elements are used) as shown in Fig. 4.7. 

The first two natural frequencies for various angles of rotation of the material axes, in 
an anisotropic cantilever beam, using both the displacement and hybrid approximations are 
tabulated in Table 4.2. The material properties and the material model assumed are the 
same as for the static anisotropic case, i.e., 3 independent constants in a crystal with cubic 
syngony, where 


E l = E 2 = 1.9716 x 10 7 lb / in 2 

^12 = 0.2875 (4.32) 

G u = 5.4758 x 10 6 lb/in 2 
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Figure 4.5: Normalized natural frequencies (mode 1 analytical solution vs. mode 1 finite 
element solution) for various quadratic finite element mesh sizes. 
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Figure 4.6: Normalized natural frequencies (mode 2 analytical solution vs. mode 2 finite 
element solution) for various quadratic finite element mesh sizes. 
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Figure 4.7: Quadratic quadralateral finite element meshes used to study the behavior ot 
hybrid stress elements in vibration analysis. 
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Table 4.2 

Natural frequencies for an anisotropic cantilever beam 
using 10 quadratic elements for different material axes orientations. 


Orientation of the Axes 

u>\ (hybrid) 

u>x (displ.) 

uj'i (hybrid) 

u> 1 (displ.) 

0° 

14.130 Hz 

14.197 Hz 

83.675 Hz 

84.682 Hz 

O 

O 

CO 

13.020 Hz 

14.061 Hz 

78.264 Hz 

83.995 Hz 

45° 

12.679 Hz 

13.058 Hz 

76.499 Hz 

78.911 Hz 

60° 

13.020 Hz 

13.707 Hz 

78.264 Hz 

82.211 Hz 

90° 

14.130 Hz 

13.276 Hz 

83.675 Hz 

80.023 Hz 


The results of the hybrid model are as expected, with frequencies falling as the angle ot 
rotation is increased, reaching a minimum at a rotation of 45°, and then increasing symmet- 
rically (since E\ = E 2 ) np to a rotation of 90°. 

The results of the displacement model do not show symmetry about 0 = 45°, and are 
not as invariant under a rotation of the axes. 

A Specific Numerical Example 

A tapered cantilever beam consisting of three crystals of the same material but with differ- 
ent orientations of the material axes is considered next and analyzed tor its displacements, 
stresses, natural frequencies, and mode shapes. The beam is shown in Fig. 4.8. Since the 
nickel alloy for which experimental data was provided exhibits cubic syngony in its crystals, 
the same material properties are considered as used in the previous. Also, since of all the 
elements tested, the 8 noded hybrid-stress element gave the best results, this element is used 
in the mesh shown in Fig. 4.9. 

The following are specifications for the beam: 

£ = 5 in; depth at fixed end = —5 in; depth at free end = 1 in; 
static load = 150 lb; crystal orientation = 0°, 30° and 45° (see Fig. 4.7). 

The results are tabulated below: 
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Table 4.3 

Comparison of hybrid displacement mode results for a 
tapered, 3-crystal anisotropic cantilever beam 

Parameter Hybrid Method Displacement Method 

Utip 0.0539 in 0.0502 in 

^bending (max) 1-632 X 10 5 psi 1.519 X lO -1 psi 

W 1.643 X 10 5 psi 1.511 x 10“ psi 

30.198 Hz 30.869 Hz 

131.482 Hz 143.976 Hz 

quencies by even less. Howevef Tthe SO^rotat ' ^ 6 " 5% ' and the natural fre “ 

mcrease rapidly to a maximum of almost 12% * t0 a 6 °° rotation ’ a11 et ' p ors 

* - « accurately 

observed in elements 1,6 while the maximum 1 f ff ^ he maximum bending stress is 

“ - - - "* 
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5 Three-Dimensional Element Definition 

to equilibrium conditions, is given by. 


7T 


me 


=£5 


r S<riV-l <T T (Du)dV + f u T Tds\ 

JV n Js<Tn ' 


(5.1) 


The stresses are interposed using the stress parameters 0 and the polynomial stress 

shape function P, (5.2) 

a — rp 

such that the homogeneous equilibrium conditions are satisfied. 

D T a = 0 


5.3) 


A primary difference between three-dimensional anc l 

calculations is in the calculation of the P matnx m the above equa • • 'functions 


n.0 < up — Mr 


5.4) 


a 1 1 VI- TO = fin and Tip = 6. so that the minimum number of 

For a quadratic twenty node bn , J ' ™ ^ be at least this minimum number of 

stress parameters is 60 - b - o4. ° - f„ n -tinns must be formed from a complete 

stress parameters, but also the stress rnterpolatron tunc ^tYqu adratics would yield 
set of basrs funettons, in this o s^X^iuhtum, the number of 

a total of 60 stress parameters but if the st ‘ e = ses ired for invariance. Thus, lull 

*• - — 

would be: 


0 


0 


P = 


0 

0 

0 

0 

0 


1 . 


0 


0 


0 


0 


0 


0 


0 

0 

0 

0 

l...z 3 

0 


0 


o.o ; 


1 ... 
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The number of stress parameters can be reduced by applying the homogeneous equilib 
rium equations: 

dcr x foxy 9 t xz 

dx ' dy dz 

dTxy dOy ^Tyz 

dx dy dz 

dr IZ dr^_ da^ 

dx dy dz 

Each of these equations yields a reduced order polynomial with 10 terms. For the equa- 
tions to be satisfied for arbitrary values of x,y, and 2 , the coefficients tor each term must 
equal zero. Thus, 30 equations are generated relating the stress parameters to each other 
and the P matrix is reduced to a 6 x 90 matrix. 


= 0 

= 0 (5.6) 

= 0 


For an isotropic case, the number of stress parameters can be further reduced bv applying 
the Beltrami-Michel stress compatibility equations. These equations are essentially a refor- 
mation of the compatibility conditions in terms of the stresses. For the anisotropic case, one 
must again start with the strain compatibility conditions and reformulate the stress compat- 
ibility equations using an anisotropic stress-strain law. The strain compatibility equations 


are: 

d 2 e x d 2 s y _ d 2 7xy 
dy 2 dx 2 dxdy 

d^ y &t a d^yz 
dz 2 + dy 2 dydz 

0 2 e z d 2 c x _ d 2 lxz 
dx 2 ^ dz 2 dxdz 


9 d 2 e x _ 

s ±( 

dlyz 

+ 

dy 

~ dydz 

dx\ 

dx 

dy 

dz 

, d% 

d ( 

&lyz 


, dy'xy 

w dxdz 

" dy \ 

dx 

dy 

T dz 

0 d 2 z z 
~ dxdy 

II 

dlyz 

dx 

r> (xz 

dy 

dz 


(5.7) 


The stress-strain law for a fully anisotropic material is given by: 

£ t = dij&j 


(5.8) 


where the compliance matrix S is 


a n 

0-12 

Ol3 

Ol4 

Ol5 

0-12 

O 22 

O 23 

024 

O 25 

0 13 

O 23 

a 33 

034 

035 

a 14 

O 24 

.0 34 

O 44 

045 

015 

025 

O 35 

O 45 

055 

aie 

026 

036 

O 46 

056 


Ol6 

0-26 

a 36 

since a t} = etj,- (5.9) 

U 46 

a 56 
0-66 
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£i = SijPjkflk 


The stress compatibility equations become 

d^SPy^k d 2 SP 2k 0k _ ^SPy^ 

dy 2 dx 2 dxdy 


(5.10) 


Since S is constant, the compliance matrix can be pulled out ot the differential, giving 

r , d 2 P lk h , „ d 2 P 2k 0 k _ c d 2 P 4k h 

^ ~~Etf~ + b ' 2 dx 2 dxdy f51] 


This procedure yields six equations, each with four terms. Again, since these equations 
must be satisfied at all points in the element, each coefficient must be zero and 24 equations 
are generated to eliminate stress parameters, leaving a total of 66 stress paiameteis. 

The solution for the P matrix can be obtained in closed form for both the constraints 
due to equilibrium and due to stress compatibility, however, the algebra is quite tedious. 
This procedure can be accomplished by a series of matrix manipulations on the original P 
matrix and on the various derivatives of the P matrix. Both the first and second ordei 
derivatives must be defined beforehand, and since the position ot a particular term (i.e., the 
xy 2 term) will vary with which derivative is being taken, a careful account of the variables 
associated with each term in the derivatives is necessary so that the constraint equations 
can be properly formulated. The P matrix is statically condensed using first the equilibiium 
constraints and then the stress compatibility constraints resulting in a 6 x 66 matrix. 

Because the number of stress parameters is above the minimum number of 54, the solution 
may be overly stiff, but this can be determined by a comparison ot this fomulation with the 
54 stress parameter element of Rubenstein, Arluri, and Punch [20] while using an isotropic 
material model. 

For an eight node brick, rip = 24 while tir = 6, so that the minimum number of stress 
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parameters up is 18. If trilinear stress interpolation functions are used, a is given by 

cr x = fix + fi2 x + #32/ + Paz 
&y — #5 + Pe X + fiiy + fis z 

& z — @9 + PiqX + fin y + fii 2 z 

(5.12) 

T xy = #13 + fi\4 x + fi\ 53 / + fiieZ 

r xs ~ Pi 7 + #18^ + #19j/ + #20- 
r y2 = #21 + fi22 x + #23 J/ + #24 2 

Introducing cr into the equilibrium equations gives 

#20 — #2 #15 

#24 = — #14 — #7 (5.13) 

#23 — #12 — #18 

Substituting and renumbering yields 

O' 'x — fix + #2^ + #32/ + #4- 

“ #5 4" #6 ^ + #72/ 4" #S- 

^3 — #9 + #10^ + #11 2/ + #12^ 

(5.14) 

r xy — #13 + #14'£ + #15?/ + #16^ 

T xz = #17 + fi\S x + #192/ ” (#2 + #15)2 

x yz — #20 + #21^ “ (#12 + fixs)y (#7 + fil4) z 


Twenty-one stress parameters remain and thus the element should be free from spurious 
zero-energy modes. 

A minimum stress parameter element from an eight noded brick has been proposed by 
Punch and Atluri [18]. Here certain symmetries of the elements are assumed and methods 
of group theory are used to eliminate some of the stress parameters without affecting the 
rotational invariance of the element. This method is clearly adequate for isotropic materials 
and may be adequate for materials which, while not isotropic, have high degrees of symmetry, 
but they are not applicable to fully anisotropic materials. 
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The stress approximation for this 18/? element is given by: 

g x — (A + A) + 2 fax + /3isi/2r 

&y = (A + A A) + 2/?8j/ + (3\jXZ 

g z — (A — A) + 2{3qz + ,3\exy 

(5.15) 

T xy = Ao + (/?14 “ A) x + (/?13 — )2/ + (A 3 + A) 2 

r x2 — Ai + (As “ A) ^ + (A + A)y ~ (A3 + AA 

r yz = A2 + (A - A - A A - (A + A)y - ( A + As A 

The stiffness matrix is calculated as before 

K = G T H~ l G (5.16) 

where 

H = [ 1 / [ P T SP\J\d£dT)d<; (5.17) 

Jo Jo Jo 

and 

G = /' I' f P T B” d^drjds (5.18) 

Jo Jo Jo 

The calculation of the compliance matrix is made in the material axes and a rotated 
compliance matrix S is calculated with respect to the global axes. The new terms are given 

b y ; 

Ojj = ^ ^ ' Q'mnQimQjn (5.13) 

m n 

where a mn are the components of the compliance matrix in the material axis and the q im 
and q jn are functions of the direction cosines relating to the material and global axes. 
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6 Alternate Hybrid Stress Element Formulations 

6,1 The Eight— Node Punch and Alturi Brick Element 

An alternative approach to the assumed stress field was proposed by Punch and Alturi. The 
basic hybrid stress model is developed in the same manner as previously outlined regardless 
of the form of stress interpolation. The Punch and Alturi approach is based upon symmetry 
group theory. The method defines a set of ‘natural’ irreducible strain subspaces which are 
invariant to element translation and to the 24 symmetric rotations of a cube. For each such 
strain space, at least one stress space must be defined which is also invariant. For an eight 
node brick, there are eighteen natural strain subspaces along with six rigid body modes. A 
minimum of 18 independent stress subspaces are required. If complete quadratic functions 
are used to generate the equilibrated stress subspaces, a total of 48 stress subspaces are 
created. Using all 48 subspaces would be equivalent to using complete polynomials and 
eliminating 12 of the resulting 60 stress parameters by applying equilibrium constraints. 

The difference with this formulation and those given earlier is that here each stress param- 
eter is related to an invariant subspace rather than to an individual term in an equilibrated 
complete polynomial. Thus, any stress parameter can be removed without affecting the in- 
variance of the resulting element. After the parameter is deleted, the rank of the resulting 
stress function matrix, G, must be at least 18. Certain linear terms must be included, but 
most of the quadratic terms can be eliminated without affecting the rank of G. However, 
once subspaces are removed, the resulting G matrix is composed of incomplete polynomials, 
and may no longer contain all of the cardinal stress states of pure bending. Therefore, the 
resulting element must be evaluated to be sure that the bending stiffness is adequate. 
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For a least order formulation, these stress subspaces are presented below as second order 
tensors in 3-space: 
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The corresponding stress space function matrix is given below: 
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This element performs well in tension, pure shear, and pure bending with isotropic ma- 
terial properties. It also performs well in both tension and pure shear with fully anisotropic 
material properties as long as the loads are applied separately and along the axis of a cubic 



element. When the loads are combined, the performance begins to degrade. However, exact 
solutions of the stresses and strains in these combined loading situations must be completed 
before the full extent of the degradation can be determined. 

One additional problem exists with this element. The formulation of the stress and strain 
subspaces was based upon the symmetric transformations of a cube. As the geometry of the 
element is distorted, the performance of the element is degraded. A similar degradation may 
occur if the anisotropy of the material is increased, especially when non-svmmetric loads are 
applied. It should be noted, though, that in most cases these elements give better results 
than eight noded displacement elements, although the computational effort is significantly 
greater. 


6.2 The Twenty Node Punch and Alturi Brick Element 

A twenty node hybrid element has been developed based upon the work of Punch and Atluri. 
This element uses the same lines and quadratic terms that are used for the eight node element 
described above, but adds six linear terms, eighteen quadratic terms, and twelve cubic terms 
for a total of 54 stress parameters. The following stress subspaces make up the twenty node 
elements: 
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The twenty node hybrid element 
The /f matrix is given by 


has an H matrix of rank 54 which must be invertible. 

H = pT S p 
1=1 


where L is the number of integration points. The matrix S is of rank six so that a minimum 
of nine integration points is required for H to be of full rank and thus invertible. Without 
a priori knowledge of the mesh, the integration scheme must be symmetric. Using a Gaus- 
sian quadrature, this requires 27 integration points (i.e., 3 x 3 x 3), and consequently the 
computational effort is 3 times what is actually required. Irons has implemented on a 14 
point integration scheme, which is basically a 2 x 2 x 2 scheme with an additional integration 
point at the center of each face of the hexahedral element. This integration scheme gave 
identical results to five decimal places to the Gaussian quadrature scheme with essentially 
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60 % of the computational effort. 

“ — 

54 term model is shown below. 
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6.3 The 42 Parameter Hybrid Stress Brick Element 

An alternative element which does not suffer from these symmetry-related problems was 
also developed. Complete polynomial stress functions are used and equilibrium constraints 
are applied, yielding a 48— parameter element. The strain compatibility equations are then 
satisfied using a fully anisotropic material model. This constraint reduces the number of 
stress parameters to 42. This is an inordinately large number of stress parameters and may 
yield an overly stiff element, but it may also yield an element that is far less sensitive to 
distortion and anisotropy. 

While the calculations for the polynomials can be derived in closed form, the algebra 
is extremely tedious. However, the constraint equations can be formed as a set of matrices 
showing the initial relationship between the stress parameters as derived from the equilibrium 
equations and the stress compatibility. 

The equilibrium constraint matrix contains only constants while the compatibility con- 
straint matrix contains ratios of the compliance matrix constants. The constraint matiices 
are first internally reduced eliminating the constrained parameters from the constraint ma- 
trices. The P matrix initially contains only the complete basis functions as: 


P 0 0 

OP 0 
0 0 0 

( 6 . 6 ) 

0 0 ... 0 

0 0 0 

0 0 P 

where P represents the following 10 terms for the eight node element: 

1 2x 2 y 2 z x 2 2 xy 2 xz y 2 2 yz z 2 (6-") 

and each 0 represents 10 zeros. For the eight node element, P is initially a 6 x 60 matrix, 

while for the twenty node element, full cubic basis functions are required so that the P 

matrix is 6 x 120. For the eight node element, there are 12 equilibrium constraint equations 
and six compatability constraint equations, while for the twenty node element there are 30 
equilibrium constraints and 24 compatibility constraints. 

The elements of the P matrix are functions of position, therefore, the P matrix must be 
calculated for each integration point and then the constraint equations are used to eliminate 
the constrained stress parameters. The P matrix is then statically condensed to yield a 
6 x 42 matrix for the eight node element and a 6 x 66 matrix for the twenty node element. 
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6.4 Numerical Experiments 

The eight node and twenty node hybrid elements were compared to the standard displace- 
ment elements for both single elements and for a six element beam. The single elements were 
tested in pure tension, pure shear, bending, and torsion using both isotropic and anisotropic 
properties. As is shown in Table 6.1, all elements gave exact solutions for pure tension and 
pure shear with isotropic material properties. For pure bending, the eight node displacement 
element is overly stiff but all of the hybrid elements again gave exact results. None of the 
elements is able to give exact results for torsion, but the hybrid elements perform as well or 
better than the corresponding displacement element. 

Table 6.1. Displacements Produced by Cardinal Stress States for Isotropic 
Material Properties 


Element 

Pure 

Tension 

Pure 

Shear 

Pure 

Bending 

Pure 

Tortion 

DM 8 

100 

100 

67 

84 

H 8-42 

100 

100 

100 

84 

H 8-18 

100 

100 

100 

84 

DM 20 

100 

100 

100 

95 

H 2 0-54 

100 

100 

100 

102 


As the degree of anisotropy is increased, the performance of the displacement elements 
decreases. Table 6.2 shows information equivalent to Table 6.1 for anisotropic material 
properties. Here the ratio of Eu/E , 33 is 3 and the material axis is rotated with respect to 
the element axis by the direction cosines given below: 


Direction Cosines: 


.5774 

.574 

.5774 

.7071 

-.7071 

.0000 

-.4082 

-.4082 

.8165 


Both of the displacement elements show deterioration in bending and torsion as the 
degree of anisotropy increases although the degradation is less for the twenty node element. 
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Table 6.2. Displacements Produced by Cardinal Stress States for Anisotropic 
Material Properties 


Element 

Pure 

Tension 

Pure 

Shear 

Pure 

Bending 

Pure 

Torsion 

DM 8 

100 

100 

46 

76 

H8-42 

100 

100 

100 

84 

H8-18 

100 

100 

100 

84 

DM 20 

100 

100 

97 

92 

H20-54 

100 

100 

100 

95 


These calculations were performed using double precision for all real variable calculations. 
The result for the twenty node elements were compared for three integration rules: a 4 x4 x4, 
a 3 x 3 x 3, and the 14 point rule proposed by Irons [7]. The differences in results were less 
than one percent. Spilker [22] stated that the 14 point rule produced some ill conditioning 
of H but no such ill conditioning was detected in these runs. Consequently the 14 point 
rule was used for all subsequent calculations except for occasional checks to assure that the 
results were indeed the same for the 14 point and 3x3x3 rules. 

A six element cantilever beam was analyzed using both eight node and twenty node 
bricks. These beams are shown in Fig. 6.1. The beams were analyzed both with a pure 
moment loading and a uniform end shear and tor isotropic material properties as well as a 
series of anisotropic material properties. Figure 6.2 shows the normalized tip displacement 
for a pure moment load with eight node bricks as a tunction of the degree of anisotropy 
where the degree of anisotropy is given by the ratio of the Young's moduli in the primary 
material axes. Figure 6.3 shows the normalized tip displacement of the cantilever beam for 
uniform end shear as a function of the degree of anisotropy. The hybrid stress elements are 
clearly less sensitive to the degree of anisotropy than the displacement. This is true even 
for the least-order formulations of Punch and Atluri which depend upon element symmetry 
for their formulation. Figure 6.4 shows a comparison of a z in the cantilever beam for a pure 
moment load on the end of the beam as a function of the degree of anisotropy. The stresses 
in the beam should not change as the material properties change and all stresses except cr z 
should be zero. This is clearly not the case for the displacement elements, even tor the twenty 
node brick. The eight node element actually gives better results than the twenty node brick 
because the stresses were interpolated at the 2x2x2 Gauss points which are the optimum 
points. 

The hybrid stress elements clearly give better displacements and stresses for highly 
anisotropic material properties than their corresponding displacement elements but at some 
calculational expense. The calculation times are shown in Table 6.3 for calculation of the 
six element cantilever beam. The hybrid stress elements require up to three times as long 
for the calculations but the displacement elements require at least twice as many elements 
to obtain the same degree of accuracy in o z . The accuracy in the shear stresses and in a x 
and a y , though, are still better in the hybrid stress elements. 
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igure 6.1: Three-dimensional finite element meshes of a cantilever beam 

A. Eight node brick elements 

B. Twenty node brick elements 


1.1 



Figure fi.2: Normalized tip displacement for a cantilever beam with pure moment loading 
for various three-dimensional hybrid stress elements vs. the analytical solution as a function 
of the degree of anisotropy. 


1-1 



Figure 6.3: Normalized tip displacement for a cantilever beam with uniform end shear for 
various three-dimensional hybrid stress elements vs. the analytical solution as a function of 
the degree of anisotropy. 



Figure 6.4: Normalized Bending stress for a cantilever beam with pure moment loading tor 
various three-dimensional hybrid stress elements vs. the analytical solution as a tunction of 
the degree of anisotropy. 
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Calculation Time for a Six Element Cantilever Beam 

DM8 DM20 H8-18 H8-42 

Time(sec) 67 552 108 562 


H20-54 

1422 


6.5 Conclusions 

The hybrid stress elements printed Inhere | materials in areas of high 

calculation of both displacements and stres g 1 provides increased accuracy 

stress gradients. The twenty node ^mately a three to one increase 

er the twenty node displacemen e e , , it q 10 provides much improved results 

computational time. The eight node the computational 

„.«r the standard eight node displacement element -th less ^ ^ e , ement provides 

time. The most surprising lesu , °wever, displacement element at one-fifth ol 
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the twenty node displacement element. 
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7 Numerical Examples 
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rials was implemented for non-isotropic mate- 

functions were incorporated into the SPaTg^F^PS?^'” 53 calculation antl dis Play 
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-dimensional finite element model of a ca 


i ^ in SS - ' eignt nouc 
ntilever beam. 1 


solid elements. 



Figure 7.2: Three- 

solid elements. 


.dimensional finite element model of a ca 


ntilever beam, 40 S82 eight 


Material properties 


for Cubic Syngony 


E\ = £2= £3 


= 1.9716E + 07 psi 

5.465S£ + 06 psi 
0.2875 


Table 7.1 


Anisotropic beam modeled with 10 S82 three-dimensional elements compared with 
two-dimensional solution for ten element model. 

3-D 2-D 


Material Axis Rot. 

u-tip 

max S 

u-tip 

max S 

(deg) 

(in) 

(psi) 

(in) 

(psi) 

0 

.0509 

14201 

.0512 

11532 

30 

.0606 

14242 

.0605 

12272 

45 

.0639 

14256 

.0638 

12566 

60 

.0606 

14230 

.0605 

12273 

90 

.0509 

14201 

.0512 

11532 


Table 7.2 

Anisotropic beam modeled with 40 S82 three-dimensional elements compared with 
two-dimensional solution for 40 element model. 

3-D 2-D 


Material Axis Rot. 

u-tip 

max S 

u-tip 

max S 

(deg) 

(in) 

(psi) 

(in) 

(psi) 

0 

.0508 

14618 

.0496 

11117 

30 

.0612 

14172 

.0579 

10579 

45 

.0660 

15341 

.0619 

11023 

60 

.0613 

15546 

.0579 

10579 

90 

.0508 

14618 

.0496 

11117 


62 



Table 7.3 


Anisotropic beam modeled with 40 S82 three-dimensional elements compared with 
three-dimensional solution for ten element model. 

3-D (40 elem) 3-D (10 elem) 


Material Axis Rot. 

u-tip 

max S 

u-tip 

max S 

(deg) 

(in) 

(psi) 

(in) 

(psi) 

0 

.0508 

14618 

.0509 

14201 

30 

.0612 

14172 

.0606 

14242 

45 

.0660 

15341 

.0639 

14256 

60 

.0613 

15546 

.0606 

14230 

90 

.0508 

14618 

.0509 

14201 


Table 7.4 

Anisotropic beam modeled with 40 S82 three-dimensional elements compared with 
three-dimensional solution for ten element model. 


3-D (40 elem) 3-D (10 elem) 


Material Axis Rot. 

Frequency 

Frequency 

(deg) 

(Hz) 

(Hz) 

0 

260.53 

260.69 

30 

238.24 

239.19 

45 

230.42 

233.27 

60 . 

238.48 

239.36 

90 

260.53 

260.69 
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Table 7.5 


Effect of orientation of material axes on frequencies for SSME turbine blade with material 
of cubic syngony to simulate the single crystal nickel alloy. 


Rotation about X-axis 

fl (Hz) 

f2 (Hz) 

f3 (Hz) 

(deg) 

0 

4379 

9388 

11818 

30 

4348 

9375 

11718 

45 

4334 

9318 

11676 

60 

4337 

9298 

11707 

90 

4376 

9388 

11816 

Rotation about Y-axis 

fl (Hz) 

f2 (Hz) 

f.3 (Hz) 

(deg) 

0 

4379 

9388 

11818 

30 

4378 

9918 

11820 

45 

4361 

10216 

11826 

60 

4362 

10101 

11834 

90 

4378 

9388 

11816 

Rotation about Z-axis 

fl (Hz) 

f2 (Hz) 

f3 (Hz) 

(deg) 

0 

4379 

93S8 

11818 

30 

4379 

9337 

11814 

45 

4380 

9306 

11815 

60 

4382 

9347 

11813 

90 

4380 

9393 

11816 


Material properties for Cubic Syngony 


El = E2 = E3 = 1.9716E + 07 psi 
G= 5.465SE + 06 psi 

v = 0.2875 
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first and third modes changed very little with material orientation. The largest change in the 
first mode frequency occurred with a material rotation of 45 degrees about the global x axis. 
This orientation resulted in a first mode frequency reduction of 1.0%. Material axis rotations 
about the global 2 axis produced very little change in frequency, the largest change being a 
0.9% reduction in the second mode frequency for a material axis rotation of 45 degrees. 

The fourth example consists of the previous blade example modified to incorporate the 
S82 solid elements in the base region of the blade, including the “fir tree” portion. A plot 
of the modified blade model is shown in Fig. 7.4. The same series of runs made for the 
previous configuration was made for this model. The material axes were again rotated about 
each of the principal axes and frequencies calculated for each orientation. Since the S82 
elements were used throughout the blade, the material orientations were effective from the 
tip to the base of the blade. The frequency results are summarized in Table 7.6. Once again, 
rotations of the material axes about the global y axis produced the greatest effect, with 
an 8.4% increase in the second mode frequency for a material axis rotation of 45 degrees. 
However, unlike the previous configuration, rotations about the global x axis also produced 
a significant effect. A reduction to the first mode frequency of -/.6% occurred for a material 
axis rotation of 45 degrees about x. Material axis rotations about the global 2 axis again 
produced very little change in frequency. The largest change was a 0.8% reduction in the 
second mode frequency for a material axis rotation of 45 degrees. The largest change in the 
third mode frequency was a 2.8% reduction which occurred for a 30 degree rotation about 
the a - axis. 
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Figure 7.4: Three-dimensional model of a SSME turbine blade modified to incorporate the 
S82 solid elements in the base region of the blade. 
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Table 7.6 


Effect of orientation of material axes on frequencies for SSME turbine blade with material 
of cubic syngony to simulate the single crystal nickel alloy. 


Rotation about X-axis 

fl (Hz) 

f2 (Hz) 

f3 (Hz) 

(deg) 

0 

4969 

8764 

11849 

30 

4679 

8755 

11516 

45 

4593 

8677 

11561 

60 

4678 

8644 

11754 

90 

4967 

8764 

11847 

Rotation about Y-axis 

fl (Hz) 

f2 (Hz) 

f3 (Hz) 

(deg) 

0 

4969 

8764 

11849 

30 

4S28 

9218 

119S0 

45 

4737 

9499 

12059 

60 

4769 

9426 

12045 

90 

4969 

8764 

11847 

Rotation about Z-axis 

fl (Hz) 

f2 (Hz) 

f3 (Hz) 

(deg) 

0 

4969 

8764 

11849 

30 

4957 

8722 

11837 

45 

4957 

8694 

11837 

60 

4963 

8728 

11836 

90 

4971 

8768 

11847 


Material properties for Cubic Syngony 


E1 = E2 = E3= 1.9716E + 07 psi 
G = 5.4758.E + 06 psi 

v = ' 0.2875 
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8 Calculation of Material Constants and Stress Mea- 
surements for Anisotropic Materials 

8.1 Introduction 

Anisotropic materials, in particular single crystal 

andmore applications in the aeros pace — 

The anisotropic propert J 

Th" “ elirr LoCpic materials ^nay he 

in the case of strong anisotropy. It was shown in ma- 

finite elements do not perform well m the static an * applications, 

tenals and that a hybrid finite element formulatior i is in the 

Similarly, strong anisotropy inherent in sing e cry , it can be observed that the 

design and interpretation of experimental results. For exa P ’ powers of 

calculation of material constants for anisotropic mate ^ plated values 

the direction cosines of the the ciystallograpnc ax . m i sa lia n ments of the crvstallo- 

of the material constants will be strongly dependent on ^ design of robust 

S parametric evaluate of —pic 

materials. 

The particular experiments considered here include: 

. the calculation of elastic material constants from tensile test experiments at an arbi- 
trarv orientation of the material axes and of the strain gauges 

. the calculation of stresses from strain measurements at an arbitrary onentation of the 
material axes and of the strain gauges 

arbU^^ 

"a’ T^econchasions^fthis'sensitivityTtudy are the basis for crimination of the 

co t::;i m<r: ^ 

developed. 
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8.2 Stress, Strain and Compliance for Anisotropic Materials 

8.2.1 Definitions 

In the general theory of elasticity the stress and strain measures are defined as second order 
tensors: <7 and e, respectively. The general representation of these tensors in any Carte 
coordinate system {x t }, * = 1,2,3 with base vectors e, is of the form (see reference [ ]). 


= (Jij e, <g> e 3 

(8.1) 

= e t j e, ® Bj 

(8.2) 

defined by: 


S = — 
dcr 

(8.3) 

e, ® e } ) © (e fc ®e,) 

(8.4) 


and has a representation 


8.2.2 Transformation Under Rotation of a Coordinate System 

Consider two different coordinate systems: {*} with base e, and {*,} with base Then 
the stress tensor a has two different representatrons m each of these systems, related by the 

transformation law: __ (8 51 

o ij = r xk r 3 l a k i 1 " ' 

where the elements of the rotation matrix rv, are defined by 

r -e e ( S6 > 

Using matrix notation, the above can be expressed as 

M = WPf 

where the rotation matrix [R] has the property 

[R ]- 1 = [*l r 

The components of the strain tensor are transformed according to a similar formula 

[el = [RM[R T } (8 ’ 9) 

Note that matrix notation is used here to emphasize the fact that tensors o* and e (defined 
as linear operators) remain unchanged and that only their representations change. 


(8.7; 


(8.81 
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8.2.3 Stress, Strain and Compliance - Technical Notation 


In practical applications a slightly different notation is usually introduced to represent stress, 
strain and compliance. Instead of the two-dimensional matrices, stress and strain vectors 
are introduced: 


<T = {<T,} = {<7n, 022, 033 i &13, 023/ (b.10) 

£ = {£,} = {£ll,£22 ? £ 33) £ 12) £ 13i£23} :r (8'H) 


The compliance matrix is defined by 


S = 


de 

far 


and has, in the most general case, 21 independent components. For cubic systems, which 
are of primary interest here, the number of independent constants is equal to 3. 


It should be noted that, with the technical notation introduced in the previous section, the 
stress vector, strain vector and compliance matrix are not elements of tensoi spaces defined 
on the three dimensional Euclidean space anymore. Therefore, they do not necessarilv 
obey the rules relating to objects of these spaces (tensor laws), in particular, the rules 
of transformation under rotation equations (8.7) and (8.9). Thus, in technical notation, 
boldface symbols represent vectors and matrices, not tensors. 

It can be shown, however, that the transformation of the components of the stress and 
strain vectors under the rotation of a coordinate system are represented by 


<y — Qa 


(8.12) 


and 

c = Qe (8-13) 

where the matrix Q consists of the appropriate combinations of products of elements of the 
rotation matrix R. Note that, since tensor laws are not valid for the technical notation, the 

inverse matrix Q ^ is not equal to and has to be reconstructed iiom elements of the 
T 

matrix R 1 . 

For the compliance matrix, it can be shown that it transforms according to the formula. 

S = QSQ~ l (814) 

Note that the compliance matrices S and S are non-symetric. The symmetric compliance 
matrices are obtained if new measures of strains 7 q = 2s, } are introduced into vector e. 
However, since this new strain vector requires a different transformation matrix than Q in 
(8.13), it will not be used in the computations. If needed, the symmetric compliance matrix 
can be easily reconstructed from the matrix S by the use of the formula 

s S, } = Sij i= 1,2,3, j = l ,...,6 (815) 

s S {j = 2 Sij 1 = 4,5,6 j = 1 6 
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8 2 4 Compliance Matrix in the Local Coordinate System 

compliance matrices may have as many as 21 mdependen te J^ whic h 
or as low as two independent constants (for isotropic materials). For cubic sys ^ 

are of primary interest here, there are three independent constants a lt a 2 , 

S in the following locations: 


S = 


a i 

a 2 

a 2 

0 

0 

0 

0>2 

a i 

a 2 

0 

0 

0 

a 2 

a 2 

a a 

0 

0 

0 

0 

0 

0 

\a 3 

0 

0 

0 

0 

0 

0 

2 a 3 

0 

0 

0 

0 

0 

0 

to te- 

a 

CO 


(8.16) 


so that * _ (8.17) 

S = La v 

with non-zero entries corresponding to actual locatio Tvnical forms of this matrix 

can easilv be reconstructed from the structure of the ma n. ■ P f prence fgl They 

for various crystal classes are shown (in the symmetric version > S) in reference [8], Ihev 

are also briefly presented in Appendix C of this report. 


8.3 


Evaluation of Material Constants for Anisotropic Materials 

8 3.1 Basic Formulation 

;:ril ^ it; rXl 

gauge (single measurement). In general, systems {x t }, {x,} an { ^ } 

Introducing the transformation matrices. ^ 

q _ representing transformation from {x,} to {x,} 


and 


Q( a ) _ representing transformation from {x,} to {x t } 


(8.19) 
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Figure 8.1: A typical test for anisotropic materials, 
the formula for measured stress e a can be written as: 


Q(?) • QSQ - = s* (8.20) 

where qJ", represents the first row of the matrix Q {a) . Note that the index (a) in parentheses 
indicates quantities associated with measurement number a, but not necessarily measured 
in the direction of x a . In particular, is represented in the coordinate system }. 

Introducing the locator matrix L to represent the specific crystal class, this equation 
becomes 


^ -- 


Q(i) • QLaQ- 1 *™ = e c 


( 8 . 21 ) 

or in the index notation 

Q{mQnuQjn<J { n a] L ljP a 0 = c a (8.22) 

It can be noted that the unknowns in this equation are material constants a,. For a tensile 
test, the stress vector is defined as <7 — {< 7 , 0, 0, 0, 0, 0} and the above equation can be recast 
in a “normalized” version 

(8.23) 


QtiQrr.a-jlUme = £ “ 


r( Q ) 


Assuming that a total of m independent measurements are made, the system of equations 
used to determine the material constants is of the form 


or, using index notation: 


Ka= e 


^ 


(8.24) 


(8.25) 


where /? — 1 ...n (number of material constants), a = / ...m (number of measurements) 
and 

Ka0 = Q[ a 2QmiQjiL, }P (8.26) 
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Note that a is a counter for single measurements (strain readings) If several strain 
Gauges are used on the same sample, then each gauge reading is counted as one measure- 
ment providing one piece of information (one equation in (8.24)). If several differen t samples 
are used to determine the material constants for the same crystal then measurements on 
consecutive samples are added as new equations to the system (8.24), with a being an accu- 
mula«ve counter of single measurements (strain readings). Note that this case, elements 
of the matrices Q may change from sample to sample, due to different orientations of the 

crystal axes. 

8.3.2 Calculation of Material Constants From “Too Many” Experiments 

The necessary condition for the system (8.24) to have a unique solution is that the number 
of single measurements is at least equal to the number of unknown coefficients m > n. In 
practice the number of single measurements will usually be larger, so that m > n, and the 
system (8.24) will be overdetermined (provided that the equations are linearly mdependen b 
For such systems, the exact solution usually does not exist and can be only found in the 

approximate sense. 

One of the popular methods for the solution of such systems is the least squares method, 
based on the minimization of the norm of the residual of the system (&--' ). 


\Ka — el 


(8.27) 


A detailed derivation of this method can be found in reference [5]. Here, we will present only 
a final form of the square system of equations to be solved: 


[K T K)a = K T e 


(8.28) 


or, in index notation: 


(K~ fa K y 0)a i 3 — K. 


"ya 


(T 


(a) 


(8.29) 


8.4 


Evaluation of Stress Components for Anisotropic Materials 


8.4.1 Basic Formulation 

Calculation of stresses in machine elements from strain gauge measurements ; is a . very typical 
operation in structural mechanics. However, for anisotropic materia s, i P .' 

be very sensitive to various experimental errors, particularly to misalignments of the strain 

gauges. 

A typical setup for obtaining measurements is presented in Figure (8.2). The gure 
present^ section^ an anisotropic body (say, a turbine blade) with the c directions o f the 
crystallographic axes defined by the coordinate system {x t }. At a given point on the 
of the body, a local Cartesian coordinate system {x,} is defined, with axes x x and x 2 
to the surLe and axis x 3 normal to this surface. In order to measure strains, a number 
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of strain gauges is located on the surface of the body (usually strain gauges are assembled 
into strain rosettes). In Figure (8.2), only one gauge, with the associated coordinate system 

{x-°^} is indicated. 

The stress-strain relationship in the material coordinate system {x;} is of the standard 
form: 

Scr = e (8.30) 

After transformation to the surface coordinate system {x,} the above equation becomes: 

QSQ~ 1 (t = £ (8.31) 

where Q represents transformation from {x,} to {x«}. Note that, in the absence of surface 
tractions, the stress components on the surface satisfy < 7,3 = 0 , i = 1 , 2 , 3 so that the essential 
non-zero components of stress can be organized into a vector 

W = {<7 n ,a 22 ,(Tu} (8.32) 

related to the three-dimensional version by a simple permutation 

a t = P 18&3 (8.33) 

As mentioned before, in order to calculate stress, one performs strain measurements in 
various directions on the surface. Each measurement provides one piece of information (one 
equation) in the form 

■QSQ- 1 PW=e a (8.34) 

where e a is a strain measurement in the gauge a. If the total of m strain gauges is used, 
then the stresses can be calculated from the system of equations 

KW — e (8.35) 

or 

K a 0 a (3 — e ot (8.36) 

where : 

/<•«, = Q^Q™Si,Q]tP«0 ( 8 - 37 ) 

and a = 1 , . . . , m, (3 — 1, 2, 3. From the structure of a matrix K it can be concluded that the 
calculated values of the stresses are sensitive to the precise specification of the orientation 
of material axes and the strain gauges. Thus, there may exist certain configuration of the 
strain gauges that will minimize this sensitivity and produce the most robust setups. This 
question will be addressed in the next section. 
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8.5 Optimization of The Strain Gauge Orientation in Tests for 
Anisotropic Materials 

8.5.1 Problem Statement 

Consider the problem of the calculation of the material constants (8.24) or of the calculation 
of stresses from strain measurements (8.64). It can be noted, that the actual form of the 
coefficient matrix K depends on the orientation of material axes (matrix and strain 
gauges for consecutive measurements (matrices Q^). Thus, it makes sense to pose the 
following problem: 

Find a combination of coefficients (in particular, Qij and Q^) such that: 

1. The system is non-singular: 

r{K) = n (8-38) 

where r(K) is the rank of the coefficient matrix, and n is the number of unknown 
material constants. 

2. The solution of the system is the least sensitive to the variation of the orientation of 
material axes or strain gauges (in case of the calculation of stresses, only the latter is 
of interest). 

The solution to this problem is suggested by the theory sensitivity oi linear systems ol 
equations, presented, for example, in reference [5]. This theory is summarized in the next 
section. 

8.5.2 Background — Sensitivity Analysis 

Consider the system of linear equations: 

Ka = 6 (8.39) 

and suppose that the coefficient matrix K is perturbed by eF where F is an arbitrary matrix 
such that ||F|| = 1 (in the appropriate norm). Then, it can be shown [5] that the sensitivity 
of the solution a to the perturbation e satisfies the inequality: 

i|a(£) ~ a| i < k(K)p(K) + o(e 2 ) (8.40) 

l! a ll 

where a(e) denotes the solution of a perturbed system of equations and: 
k(K) = Ilk'll ||i^ -1 || = spectral radius of IT 

p K = — - — = norm of the variation of K 

\\ K \\ 
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If we decide to use the natural square norm || • || 2 , then we have ||i ^||2 — cr max (AT) and 
k{K ) = a max {K)j(r^{K) and the sensitivity inequality becomes 

H a(£ ) ~ a H 2 < I + o(e 2 ) (8-41) 

||a|| 2 - a mn (K) 

In the above, <7min(-?0 is the minimum singular value of the matrix K. 

The calculation of singular values is not a very typical operation and there is no software 
readily available for this procedure. Then, it is useful to observe that singular values of K 
are defined as square roots of the eigenvalues of K T K , namely 


a 


!*> = at * 1 


i = 1, • • ■ ,n 


(8.42) 


8.5.3 Optimization of Strain Gauge Orientation and Material Axes 

With the background presented in the previous section, it is easy to observe that in order to 
satisfy criterion (2) presented in Section 8.8l, we need to: 

Find a combination of the orientations of the strain gauges and the material axes , 
which maximizes the minimum eigenvalue of the matrix A A . 


Note, that the results of this analysis are useful for the satisfaction of criterion ( 1) - existence 

of the solution. This is because the rank of matrix K is equal to the number ol nonzero 

singular values of AT, which, in turn, is equal to the number of nonzero eigenvalues of K K 
° - - (A 1 A ) • 


is greater 


Therefore, the solution to the system exists if the minimum eigenvalue A 
than zero. 

To present the above optimization criterion in a more formal way, we denote all the 
free optimization parameters by p;, i — 1, . . . A (these may be angles of strain gauges, Millet 
indices, etc.). Then, to satisfy criteria (1) and (2), one needs to: 

Find {p t } and corresponding A^ n such that: 


-r-(I\ T K) 

''min 


max 

{p.> 


At< T K 


> 0 


(8.43) 


It should be noted that this general approach can be used both for the optimization of the 
calculation of material constants and the calculation of stresses from strains. 


8.5.4 Numerical Procedure 

In order to solve the optimization problem (8.43), a variety of procedures may be applied. 
In this work, it was assumed that the number of free parameters is small enough and the 
evaluation of eigenvalue X^ K) is cheap enough that a simple searching procedure can be 
effectively applied. This procedure is based on two steps: 
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1. Scan all the possible values of {p,} with sufficiently small resolution and evaluate, lor 
each combination, the minimum eigenvalue Aj* n A) (using, for example, the Jacobi 
method [1]). 

2. Select the combination {pj, corresponding to the maximum: 


j(K t K) 

A min 


= max 

{p.} 


(a!*:*') 


This procedure is general enough to be applied to the optimization of the orientation of 
material axes and of location and orientation of strain gauges. 


8.6 Numerical Examples 

The formulation and procedures, presented above, were the basis for the development ot a 
computer software package designed to optimize the orientation of strain gauges in tensile 
tests and in stress measurements on real samples. In this section, a few selected introductoiy 
examples of this optimization will be presented. It should be noted that the formulation 
presented in this section and implemented in the program is very general and can be applied 
to any type of crystal as well as to isotropic materials (the differences between various 
material classes are restricted to the locator matrix L). For the sake of simplicity and 
easy intuitive verification, the examples presented here deal with rather simple classes ot 
materials, namely isotropic materials and cubic systems. 


8.6.1 Optimization of the Calculation of Material Constants 

The first example illustrating the correctness of the procedure is the calculation of mateiial 
constants for an isotropic material from tensile tests (Fig. 8.3). In this case, there are 
two independent material constants, so it suffices to consider one sample with two strain 
gauges. In order to test the optimization procedure, we assume that the orientations ot both 
gauges are unknown, so there are two independent optimization parameters, namely the 
angles cq and cq>. After application of the optimization procedure, the calculated optimal 
configuration of the strain gauges is: cq = 0,a 2 = 90 (or any equivalent configuration). At 
this configuration, the calculation of material constants is the least sensitive to alignment 
errors, a result that intuitively seems to be correct. This optimal configuration is presented 
in Fig. 8.3b. The same configuration is obtained if the three gauge rosette presented in 
Fig 8 3c is considered. In this case, there is one optimization parameter (angle a) and the 
optimal orientation is a = 0 (or, equivalently, a = 90, 180, etc.). In another analysis for this 
material, a triangular strain rosette presented in Fig. 8. 3d was considered. This rosette has 
one optimization parameter, namely the rotation angle a. The results of the optimization 
analysis' indicate that if the triangle is equilateral, then all the configurations are equivalent, 
i.e., sensitivity of results to misalignments of a is the same (and rather low) at any position 
of the rosette. 
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Figure 8.4: optimization of strain gauges for cubic syngony. (a) configuration of a triangular 
rosette and (b) three independent strain gauges. 

The single crystal alloys used in the manufacturing of turbine blades are usually of cubic 
syngony. For this crystal system, there are three independent material constants, so that a 
minimum of three measurements (strain gauges) are necessary. Thus, as the first example, 
we considered one sample with the material axes aligned with the edges of a sample and 
with the direction of tension (Fig. 8.4) (Miller indices for consecutive sample axes were: 
< 100 >, < 010 >, < 001 >). The first optimization case considered was a triangular 
rosette with one optimization parameter, the angle a (Fig. 8.4a). The answer obtained 
in this case is that there exists no configuration of the rosette that can provide sufficient 
information (the system of equations is always singular). The same result was obtained even 
if three or more independent strain gauges were considered (Fig. 8.4b). The explanation 
of this result is simple: one of the independent material constants for the cubic system 
is the shear modulus and, for the perfectly aligned crystal configuration presented in Fig. 
8.4, a tensile test does not produce any shear mode of deformation. Thus, there exists no 
information available to calculate the shear modulus. 

As the next case we considered two samples of cubic structure: the one analyzed before 
and a second sample, with Miller indices corresponding to the primary, secondary and tertiary 
axes < 1,1,1 >, < 1,2,— 3 >, and < —5,4,1 >, respectively. In each sample, a two-gauge 
rosette was used, with both the location (face 1 or 2) and the orientation (angle a) considered 
as optimization parameters. The solution obtained in this case is presented in Fig. 8.5a (both 
angles are equal to zero). If three-gauge rosettes are considered, the optimal configuration is 
slightly different (see Fig. 8.5b). In this case, the optimum angles are ot\ — 20°, <22 = 0 with 
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Figure S.5: Optimization of strain gauges for cubic syngony. (a) optimal configuration of 
two-gauge rosettes and (b) optimal configuration of three-gauge rosettes. 

the location of gauges as presented in Fig. 8. 5b. Note that in order to calculate the optimal 
configuration, all the samples have to be considered simultaneously, because the results are 
"coupled.' 7 For example, if the Miller indices of sample 2 of Fig. 5b were changed, then the 
optimum orientation of the rosette in sample 1 would be different from the above example. 

These simple examples illustrate the basics of the optimization procedure and the type of 
results obtainable. As previously mentioned, it is possible with this formulation to consider 
more complex crystal classes than cubic, since the procedure developed here is completely 
general. 

8.6.2 Optimization of Calculation of Stresses from Strain Measurements 

As an example of calculation of stresses from strain measurements, let us consider a sample 
of the material of cubic syngony, presented in Fig. 8.6. 

The axes of the material coordinate system {i t } are defined by three vectors represented 


the reference coordinate system E 3 as: 


e, = {0., 

0.3162, 

-0.9487} 

e 2 = {0., 

0.9487, 

0.3162} 

e ~3 = { 1 1 

o, 

0} 

The three independent constants for 

the cubic 

material are a 2 = 1.0 x 1 0 — 6 , a 2 = —0.3 x 
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Figure 8.6: Optimization of orientation of strain rosette for the evaluation of 
anisotropic samples. 



10 -6 , a 3 = 4.0 x 10 -6 , so that the symmetric compliance matrix defined in the material 
coordinate system is of the form 


1.0 

-0.3 

-0.3 

0 

0 

0 

-0.3 

1.0 

-0.3 

0 

0 

0 

-0.3 

-0.3 

1.0 

0 

0 

0 

0 

0 

0 

4.0 

0 

0 

0 

0 

0 

0 

4.0 

0 

0 

0 

0 

0 

0 

4.0 


The coordinate system {x,}, in which stresses will be calculated, is defined by three 
vectors ei, i = 1,2,3, specified as: 


ej = {0.5774, 

-0.2573, 

0.7715} 

e2 = {0.5774, 

-0.5345, 

-0.6172} 

e 3 = {0.5744, 

0.8018, 

-0.1543} 


The strains are measured by a three gauge strain rosette, presented in Fig. 8.6. The 
orientation of this rosette is an optimization parameter in this case. 

After the execution of the optimization version of the OPTAM-S code, the optimal value 
of the angle a was calculated to be 95°, which corresponds to the orientation presented in 

Fig. 8.6. 

For the configuration obtained from the optimization run, a stress calculation version of 
the code was executed, with strain readings corresponding to consecutive gauges of the rosette 
being: t (1) = 0.01, s (2) = 0.0035, c (3) = -0.003. The calculated non-zero components of 
the stress tensor are: ffn = 296.7, <7^2 = —524.8, (T 22 = 6986.3. 
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A Appendix — Description of Spar Reference Manual 
Updates 

A.l Tab Processor Updates 

A new TAB sub-processor, SMAT, was created to provide an alternative to using AUS/TABLE 
to generate material flexibility coefficients (compliance matrices) for use with solid elements. 

Based on inputs provided by the user, SMAT generates entries in PROP BTAB 2 21 as 
required by ELD when using solid elements. SMAT may be used with the previously existing 
solid element types S41, S61, and S81 as well as with the new element types S42 and S82. 

The input required by SMAT is listed below. The detailed description of this input is 
given in Appendix B. 

n, w , nre f 
EuE 2 , E3 
G 12 •> Giz, G 23 

V\2, ^13 1 ^23 

(DIRCOS) (L J),J = 1,3) ) 

(DIRCOS) (2, J), J = 1,3) i input only if nref = 0, or blank 
(DIRCOS) (3, J), J = 1,3) , 

&x 1 ®.y ? Q-z 

Y Y Y Y Y Y 

A. 2 ELD Processor Updates 

ELD was modified to accept for the new solid element types S42 and SS2. These element 
types represent the COMCO three-dimensional anisotropic hybrid stress elements. 

The S42 and S82 element material properties are obtained from the SPAR data set PROP 
BTAB 2 21 as with the other solid element types. This data set may be constructed with 
the AUS processor as before or with the new TAB subprocessor, SMAT. 

The system diagonal mass matrix, DEM, produced by the processor E includes terms 
for S42 and S82 element types. No consistent mass matrix is available for any of the solid 
element types at this time. 

The same mesh generation capabilities exist for the S82 element as tor the S81 element, 
i.e., the hexahedral element network generator. 
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A. 3 EKS Processor Updates 

, to EKS the element intrinsic stiffness and stress matrix generator, 
Modifications were made to LJXb, me eiem 

to incorporate the S42 and S82 element types. ^ ^ of the SPA R updates 

Although the corde was modified, adapted to the SPAR format 

were associated with EKS_ T intrinsic stiffness matrices inserted into the element 

the system stress matr, generator. No 

modifications to K were required. EFIL ” and “SS2 EFIL,” which contain 

The E processor now generates the data se lements respectively. E also generates 

the element information packets for theS4 ad ^82 ms for S42 and SS2 elements, 

the svstem diagonal mass matrix, DEM, which 
However, no modifications were required to b. 


A 4 GSF/PSF Processor Updates 


, , . fnr t u e n ew S42 and S82 elements. Output 

GSF was updated to generate stress data sets for the respectively, which 

hCTRQ C 49 set neon and bito oo- 
data sets are nam produced for other element types. 

is consistent with data mav be included in GSF input as de- 
control cards designating S4 and - - of this kind are not given, 

scribed in the SPAR Reference Manual R hi ^ ^ ^ ^ elements 

"Tuples lllde to PSF. A btlTcTCeb' the 

T I^TT“m=em are now ev^t the element corners insteac 
of the integration points as was the case in the earlier 

A 5 Plot Processor Updates 

GGS and PLTTEX were updated to recognize and plot S42 and ^ebmenta These 

elements may be selected by element type, group, and index, 

specified. 

A. 6 EADS Data Sets 

The SPAR modifications have been fully implemented on the NASA/MSFC EADS computer 
Jyslem "11 source codes reside in the follow, ng data sets: 
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TAB CMLP410.SPAR.FORT(PRGTAB) 

EDS CMLP410.SPAR.FORT(PRGELD) 

EKS CMLP410.SPAR.FORT(PRGEKS) 

GSF CMLP410.SPAR.FORT(PRGGSF) 

GGS CMLP410.SPAR.FPRT(PRGGGS) 

PLTTEX CMLP410.SPAR.FORT(PLTTEX) 
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B Appendix— SPAR Reference Manual Updates 



APPENDIX B - SPAR REFERENCE MANUAL - ( CONTENTS ) 

Section Foreward 

1 Introduction 

1.1 NewUserOrientation 

1.2 SPAR Overview 

2 Basic Information 

2. 1 Reference Frame Terminology 

2.2 The Data Complex 

2.3 Card Input Rules 

2.3.1 Equivalence of Word Terminators 

2.3.2 Continuation Cards 

2.3.3 Loop-Limit Format 

2.4 Reset Controls, Core size control, and the Online command 

2.5 Data set structure 

2.5.1 TABLE 

2.5.2 SYSVEC 

2.5.3 ELDATA 

2.6 Error Messages 
3 Structure Definition 

3.1 T AB - B asic Table Inputs 

3.1.1 Text 

3.1.2 Material Constants 

3.1.3 Di stributed Weight 

3.1.4 Alternate Reference Frames 

3.1.5 Joint Locations 

3.1.6 Joint Reference Frames 

3.1.7 Beam Orientation 

3.1.8 Beam Rigid Links 

3.1.9 E21 Section Properties 

3.1.10 E22, E25 Section Properties 

3.1.11 E23 Section Properties 

3.1.12 E24 Section Properties 

3.1.13 Shell Section Properties 

3.1.14 Panel Section Properties 

3.1.15 Constraint Definition 

3.1.16 Joint Elimination Sequence 

3.1.17 Rigid Masses 


(MATC) 

(NSW) 

(ALTREF) 

(JLOC) 

(JREF) 

(MREF) 

(BRL) 

(BA) 

(BB) 

(BC) 

(BD) 

(SA) 

(SB) 
(CON) 
(JSEQ) 
(RMASS) 
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3.1.18 Solid Element Materials 


(SMAT) 


TABLE 1-1: SPAR Element Repertoire. 

See Volume 1 

NAME 

DESCRIPTION 

Sections: 

E21 

General straight beam elements such as channels, 
wide-flanges, angles, tubes, zees, etc. 

3. 1.7-9 

E22 

Beams for which the intrinsic stiffness matrix is 
given. 

3.1.10 

E23 

Bar - Axial stiffness only. 

3.1.11 

E24 

Plane beam. 

3.1.12 

E25 

Zero-length element used to elastically connect 
geometrically coincident joints. 

3.1.10 


Two-dimensional (area) elements: 

3.1.13 

E31 

Triangular membrane 


E32 

Triangular plate. 


E33 

Triangular combined membrane and bending element. 


E41 

Quadrilateral membrane. 


E42 

Quadrilateral plate. 


E43 

Quadrilateral combined membrane and bending element. 

3.1.14 

E44 

Quadrilateral shear panel. 


Three-dimensional solids: 

3. 2.2. 3 

S41 

Tetrahedron (pyramid). 


S42 

4— node tetrahedral hybrid (COMCO) 


S61 

Pentahedron (wedge). 


S81 

Hexahedron (brick). 


S82 

8-node anisotropic hybrid (COMCO) 

12., 3. 2.2.3 


Compressible fluid elements: 

F41 

Tetrahedron (pyramid). 


F61 

Pentahedron (wedge). 


F81 

Hexahedron (brick). 



Notes: 

- See Section 7.2 for examples of stress output. 

- See Volume 2 (Theory) for element formulation details. 

- Aeolotropic constitutive relations permitted, all area elements. 

- Laminated cross sections permitted for E33, E43. 

- Membrane/bending coupling permitted for E33, E43. 

- E41, E42, E43, E44 may be warped. 
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- Aeolotropic constitutive relations permitted for 3-D solids. 

- Non-structural mass permitted for line and area elements. 
3.1.18 Solid Element Materials (SMAT) 


Based on inputs provided by the user, SMAT generates entries in PROP BTAB 2 21 as required by 
ELD when using solid elements. A description of the contents of each input record to SMAT 
follows. 

n, w, nref 

E„ E 2 , E 3 

G 12’ (-*13, G 2 3 
V 12’ V 13> V 23 

(DIRCOS (1J), J - 1,3) 

(DIRCOS (2J),J = 1,3) input only if nref = 0, or blank 

(DIRCOS (3J),J= 1,3) 
a x a y a z 

Y Y Y Y Y Y 

1 xx, 1 yy, 1 zz, xy , 1 yz, 1 zx 

where 

n = the material constant entry 
w = weight density (weight/unit volume) 
nref = Alternate Reference Frame (see ALTREF) 

If nref > 0, frame nref specifies material orientation relative to the element axes. 

If nref - 0, this orientation is given by DIRCOS values. 


E, = Modulus of elasticity in material direction - 1 
E 2 = Modulus of elasticity in material direction - 2 
E 3 - Modulus of elasticity in mateial direction - 3 

G 12 = Shear modulus in 1-2 plane 
G i3 = Shear modulus in 1-3 plane 
G 23 = Shear modulus in 2-3 plane 


v 12 = Poisson’s ratio of comp in dir— 2 to tension in dir— 1 

v 13 = Poisonn’s ratio of comp in dir-3 to tension in dir-1 

v 23 = Poisson’s ratio of comp in dir— 3 to tension in dir— 2 


DIRCOS (I,J) = Direction cosine matrix relating the material axis-i to the element axis-j. 
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a x , a , a z = Linear thermal expansion coefficients 

Y Y Y Y Y Y = reference or yield stress for use in stress displays. 

1 xx 1 yy 22 1 xy x yz zx J 

See PSF. 

EXAMPLE: 


Z 



a 3 = G 23 = 15 X IQS 


V -2 = 


= .3 


SMAT: 1 .283 

30. +6, 30. +6, 30. +5 
15. +6, 15. +5, 15. +5 
0.3, 0.3, 0.3 
0.866,0.000,0.500 
0 . 000 , 1 . 000 , 0.000 
-.500,0.000,0.866 
0 ., 0 ., 0 . 

l.,l.,l.,l.,l.,l. 


or 

ALTREF: 2 2,-30. 

SMAT: 1 .283 2 

30. +6, 30. +6, 30. +5 
15.+6,15.+5,15.+5 
0.3,0. 3, 0.3 
0 ., 0 ., 0 . 

l.,l.,l.,l.,l.,l. 
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3. 2. 2. 3 Three-Dimensional Elements 


Only one table pointer, NSECT (or the synonym NPROP), applies. The default value of 
NSECT is 1. For fluid elements, NSECT points to a line in a table named PROP BTAB 2 20. 
For solids, NSECT points to a line in PROP BTAB 2 21. Before executing ELD, the user 
must construct these tables via AUS/TABLE, as indicated below. Mesh generation facilities 
are described at the end of this section. 

Fluid elements F41, F61 , F81: 

For additional information see Section 12. It should be noted that FSM is the only processor 
which produces system matrices containing fluid element terms. Fluid element terms are not 
included in the system diagonal mass matrix, DEM, produced by processor E, nor in the 
system matrices produced by K , M, or KG. No form of static temperature, dislocation^, or 
pressure loading is defined for fluid elements. GSF produces no stress data for fluid elements 
Section properties are defined as follows: 

@XQT AUS 

TABLE (TV/ = 2, NJ = the number of different fluids): PROP BTAB 2 20 
J = 1 : p, $ Mass density, bulk modulus for fluid 1. 
j -2 : p, (5 $ Mass density, bulk modulus for fluid 2. 

Solid elements S41, S42, S61, S81, S82 : 

Solid element terms are included in the system diagonal mass matrix, DEM, produced by E, 
and in the system matrices produced by K and M , but not those produced by KG . Properties 

are defined as follows: 

@XQT AUS 

TABLE (NI = 31 , NJ = number of different solids): PROP BTAB 2 21 
J - 1$ Properties of material 1 follow. 

w> 

a„> 

3>/ @2 

a 3l a 32 a 33> 

Q 41 0*2 0-43 3*4 > 

35 1 3 52 3SJ 

35/ 352 3$J 354 355 °66 > 
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a x a y a z > 

Y xx Y yy ^zz Y xy ^ yz ^ zx$ 

J =2$ Properties of material 2 follow. (Same sequence as above). 
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Appendix C 

Compliance Matrices for Various Crystal Classes 
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C Appendix — Compliance Matrices for Various Crys- 
tal Classes 

In this appendix, the specific forms of compliance matrices and independent material con- 
stants are presented. The compliance matrices presented here are of the symmetric form 
s S is obtained when the engineering strain measures A t; = ± j are used m the strain 

vector. Note that the nonsymmetric compliance matrices S, corresponding to the use of 
strain measures £,j, can be easily calculated according to the formula. 

s i3 = s Sij i = 1,2,3 j = 1 , . . . , 6 

Sij = \ s S a i = 4,5,6 j = 1 6 

The independent material constants and compliance matrices for various types of crystal 
classes are listed below: 


TRICLINIC - classes 1,2 

For the material of triclinic structure in classes 1,2, there exists 21 independent compliance 
constants a u . . . ,a 21 . The location of these constants in the symmetric compliance matrix 

is defined by: 


5 S = 


a 2 a 3 

a 4 


<*6 

a 7 a 8 


a io 

an 

a 12 

<*13 

a\ 4 

<3 15 



a it 

«18 



a is 

&20 


«21 J 


MONOCLINIC - classes 3-5 

For the material of monoclinic structure in classes 3-5, there exists 13 independent com- 
pliance constants fll , . . . ,a 13 . The location of these constants in the symmetric compliance 
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matrix is defined by: 



CL\ CLj 

a s 

0 

0 

a 9 

a 2 

a lQ 

0 

0 

a n 



0 

0 

a 12 



&4 

“13 

0 


ORTHORHOMBIC - classes 6-8 

For the material of orthorhombic structure in classes 6 - 8 , there exists eight independent 
compliance constants a \. . . . . ct s . The location of these constants in the symmetric compliance 
matrix is defined by: 


’5 = 


a i a- a-j 

a 2 cig 

a 3 


0 

0 

0 

a 4 


0 
0 
0 
0 
a 5 


0 

0 

0 

0 

0 

a 6 


TETRAGONAL - cl asses 9-11 


For the material of tetragonal structure in classes 9-11, there exists seven independent 
compliance constants , . . . , ci 7 . The location of these constants in the symmetric compliance 
matrix is defined by: 


5 S = 


a\ a { 3 0 0 a 7 

a,} a 6 0 0 — a 7 

a 2 0 0 0 

&3 0 0 

a 3 0 


U4 


TETRAGONAL - classes 12-15 

For the material of tetragonal structure in classes 12-15, there exists six independent com- 
pliance constants The location of these constants in the symmetric compliance 
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matrix is defined by: 


’ S = 


a,i a 5 a 6 0 0 

a\ a 6 0 0 

0-2 


0 0 
a 3 0 


TRIGONAL-HEXAGONAL - classes 16-17 


0 

0 

0 

0 


a 3 0 

a 4 


compliance matrix is defined by: 


; S = 


, ag- The location of these c 

CI 5 Cl$ 

-O 7 

0 

a 5 —a ,6 

as 

0 

a 0 0 

0 

0 

«3 

0 

2 o 8 


&3 

2 o 6 



2 (ai - o 4 ) 


TRIGONAL-HEXAGONAL - classes 18-20 

For the material of trigonal-hexagonal structure in classes 18-20, there exist six independent 
compliance constants a lt ...,a 6 . The location of these constants in the symmetric compliance 
matrix is defined by: 

0 
0 
0 
0 

2 06 

2 (oi - o 4 ) 

TRIGONAL-HEXAGONAL - classes 21-27 

For the material of trigonal-hexagonal structure in classes 21-27, there exist five independent 
compliance constants o a , . . . , 1 5 . The location of these constants in the symmetric compliance 


’5 = 


Oi 04 Os 06 0 

U\ Cl 5 — 06 0 

02 0 0 

a 3 0 
a 3 
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matrix is defined by: 


s S — 


a x a 4 a 5 0 0 

ai as 0 0 

a 2 0 0 

a 3 0 


a 3 


0 

0 

0 

0 

0 


2(ai - a 4 ) 


CUBIC - classes 28-32 


For the material of cubic structure there exist three independent compliance constants a x .a 2 
and a 3 . Their inderpretation in terms of Young modulus E, Poisson's ratio v and shear 
modulus G is given by: 


ai 


a 2 


1 

E 


1 

Ve 


«3 = 


1 

G 


The location of material constants in the symmetric compliance matrix is defined by 


a x et2 a 2 0 0 0 

a 4 a-2 0 0 0 

0 0 0 

a 3 0 0 

0,3 0 


03 


ISOTROPIC 

For the isotropic material there exist two independent compliance constants a x and a 2 . Their 
interpretation in terms of Young modulus E and Poisson’s ratio u is: 

1 

ai = E 

1 
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The location of material constants in the symmetric compliance matrix is defined by: 

ai a 2 a 2 0 0 0 

a.i a 2 0 0 0 

ai 0 0 0 

2(ai — a 2 ) 0 ^ 

2(ai - a 2 ) 0 

2(ai - a 2 ) 

A more detailed discussion of crystal structure and compliance matrices for various crystal 
classes can be found in reference [8]. 

LOCATOR MATRICES 

The non-zero elements of the locator matrices L for each of the compliance matrices can 
be easily reconstructed by considering consecurive independent material constants a k an 
observing that s L ijk is the coefficient in matrix s S i3 corresponding to the constant a k . tor 
example, for isotropic material, the "layer” of the compliance matrix corresponding to a, is 

1 0 0 0 0 0 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 2 0 0 

0 0 0 0 2 0 

0 0 0 0 0 2 

and the "layer” corresponding to a 2 is 

0 110 0 0 

1 0 0 0 0 0 

1 1 0 0 0 0 

0 0 0 -2 0 0 

000 0 -2 0 

0 0 0 0 0 -2 
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